home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Software Vault: The Gold Collection
/
Software Vault - The Gold Collection (American Databankers) (1993).ISO
/
cdr53
/
acmalg01.zip
/
ACM633.FOR
< prev
next >
Wrap
Text File
|
1993-01-01
|
109KB
|
3,541 lines
C ALGORITHM 633 COLLECTED ALGORITHMS FROM ACM.
C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 2,
C JUN., 1985, PP. 170-182.
INTEGER I, IPRINT, J, KRES, LDL, LDX, LMAX, LOCIW, LOCW,
* N, P, P2
INTEGER ITITLE(72), IWORK(50), LVAR(6,10)
DOUBLE PRECISION PCERR, WORK(500), X(20,6), XC(20,6)
LDX = 20
LOCW = 500
LDL = 6
LOCIW = 50
READ (5,99999) (ITITLE(I),I=1,72)
WRITE (6,99996) (ITITLE(I),I=1,72)
DO 60 K=1,2
READ (5,99998) N, P, JOBAB, PCERR, IPRINT, KRES
IF (K.GE.2) GO TO 20
DO 10 I=1,N
READ (5,99997) (X(I,J),J=1,P)
10 CONTINUE
20 CONTINUE
DO 40 I=1,N
DO 30 J=1,P
XC(I,J) = X(I,J)
30 CONTINUE
40 CONTINUE
LMAX = 1
IF (K.EQ.1) GO TO 50
P2 = 2
LMAX = 6
50 CALL LDA(XC, LDX, N, P, JOBAB, PCERR, P2, LVAR, LDL,
* LMAX, KRES, IPRINT, WORK, LOCW, IWORK, LOCIW)
60 CONTINUE
STOP
99999 FORMAT (72A1)
99998 FORMAT (3I3, F3.0, 10I3)
99997 FORMAT (F5.1, F7.0, F5.0, F5.0, F7.0, F5.0, F6.0)
99996 FORMAT (1H1, 72A1)
END
C
C --------------------------------------------------------------
C
SUBROUTINE LDA(X, LDX, N, P, JOBAB, PCERR, P2, LVAR, LDL,
* LMAX, KRES, IPRINT, WORK, LOCW, IWORK, LOCIW)
C ***
C
C LDA PERFORMS A LINEAR DEPENDENCY ANALYSIS, WHICH ATTEMPTS TO
C IDENTIFY THE NUMBER OF LINEAR DEPENDENCIES AND THE BEST SUBSETS
C OF ESTIMATED AND PREDICTOR VARIABLES.
C
C DESCRIPTION OF ARGUMENTS:
C
C X REAL(LDX,P). LDX .GE. MAX(N,P).
C MATRIX CONTAINING THE N OBSERVATIONS IN ROWS WITH THE
C P VARIABLES IN COLUMNS. X IS DESTROYED BY LDA.
C
C LDX INTEGER.
C LDX IS THE LEADING DIMENSION OF THE ARRAY X. LDX MUST
C BE GREATER THAN OR EQUAL TO THE MAXIMUM OF N AND P.
C
C N INTEGER. N .GT. 1
C N IS THE NUMBER OF ROWS (OBSERVATIONS) IN THE MATRIX X.
C
C P INTEGER. P .GT. 1
C P IS THE NUMBER OF COLUMNS (VARIABLES) IN THE MATRIX X.
C
C JOBAB INTEGER.
C JOBAB CONTROLS THE TYPE OF COMPUTATIONS PERFORMED.
C IT HAS THE DECIMAL EXPANSION AB WITH THE FOLLOWING
C MEANING:
C
C A = 1 DETERMINE NUMBER OF DEPENDENCIES FROM
C PCERR.
C = 2 NUMBER OF DEPENDENCIES IS SPECIFIED BY
C INPUT PARAMETER P2.
C = 3 NUMBER OF DEPENDENCIES IS SPECIFIED BY
C P2 AND SPECIFIED SETS OF ESTIMATED
C VARAIABLES GIVEN IN LVAR.
C
C B = 0 COMPUTE ONLY BENCHMARK SOLUTION (USUALLY
C CLOSE TO OPTIMAL SOLUTION)
C = 1 TRY ALL POSSIBLE SUBSETS OF VARIABLES
C TO FIND SETS WITH SMALLEST "REDUCED
C SPACED" RESIDUALS.
C FOR EXAMPLE: IF JOBAB = 10, LDA WILL USE PCERR TO
C DETERMINE THE NUMBER OF DEPENDENCIES AND ONLY
C COMPUTE THE BENCHMARK SOLUTION.
C
C PCERR REAL. PCERR .GE. 0 AND PCERR .LE. 100.
C WHEN A = 1 IN JOBAB, PCERR IS THE MAXIMUM PERCENTAGE
C OF UNEXPLAINED ERROR VARIATION WHICH IS ACCEPTABLE,
C ASSUMING A DIAGONAL STRUCTURE IN THE ERROR MATRIX.
C IF A = 2 OR 3 IN JOBAB, PCERR IS NOT REFERENCED.
C
C P2 INTEGER. P2 .GT. 0 AND P2 .LT. P.
C IF A = 2 OR 3 IN JOBAB, P2 IS THE NUMBER OF DEPENDENT
C VARIABLES FOR THE ANALYSIS.
C IF A = 1 IN JOBAB, P2 MUST APPEAR AS AN ARGUMENT BUT
C MAY BE UNSPECIFIED. THE COMPUTED VALUE OF P2 IS
C RETURNED.
C
C LVAR INTEGER(LDL,LMAX). LDL .GE. P2.
C IF A = 1 OR 2 IN JOBAB, LVAR MUST APPEAR AS AN ARGUMENT
C BUT MAY BE UNSPECIFIED.
C IF A = 3 IN JOBAB, LVAR CONTAINS LMAX SETS OF ESTIMATED
C VARIABLES TO ANALYZE. EACH SET CONTAINS P2
C VARIABLES.
C
C LDL INTEGER.
C LDL IS THE LEADING DIMENSION OF THE ARRAY LVAR. LDL
C MUST BE GREATER THAN OR EQUAL TO THE INPUT VALUE OF P2
C (A=2,3) OR ITS EXPECTED COMPUTED VALUE (A=1).
C
C LMAX INTEGER.
C FOR A = 1 OR 2 IN JOBAB:
C IF B = 1 IN JOBAB, LMAX IS THE MAXIMUM NUMBER OF
C COMPETING DEPENDENCIES TO OUTPUT FOR USER
C INSPECTION.
C IF B = 0 IN JOBAB, LMAX MUST BE SET TO 1.
C FOR A = 3 IN JOBAB:
C LMAX IS THE NUMBER OF SETS OF ESTIMATED VARIABLES
C IN LVAR TO ANALYZE.
C
C KRES INTEGER.
C IF B = 1 IN JOBAB, A SET OF ESTIMATED VARIABLES WILL
C NOT BE ANALYZED FURTHER IF THE "REDUCED SPACE"
C RESIDUAL IS 10**KRES TIMES LARGER THAN THE SMALLEST
C CURRENT "REDUCED SPACE" RESIDUAL COMPUTED (BENCHMARK
C SET EVALUATED FIRST). IF KRES IS LESS THAN ZERO, THE
C DEFAULT VALUE OF 2 IS USED.
C IF B = 0 IN JOBAB, KRES IS NOT REFERENCED.
C
C IPRINT INTEGER. IPRINT .GE. 1 AND IPRINT .LE. 4.
C IPRINT CONTROLS THE AMOUNT OF OUTPUT. FOR MOST ANALYSIS
C IPRINT = 1 WILL PROBABLY BE SATISFACTORY. AS IPRINT
C INCREASES, THE AMOUNT OF OUTPUT INCREASES.
C
C WORK REAL(LOCW).
C SCRATCH STORAGE USED BY LDA. USER NEED NOT SPECIFY
C ANY ELEMENTS IN WORK.
C
C LOCW INTEGER.
C THE NUMBER OF ELEMENTS IN WORK. LOCW MUST BE LARGER
C THAN N*P + 3*P**2 + 4*P + N + 5*LMAX.
C
C IWORK INTEGER(LOCIW).
C SCRATCH STORAGE USED BY LDA. USER NEED NOT SPECIFY
C ANY ELEMENTS IN IWORK.
C
C LOCIW INTEGER.
C THE NUMBER OF ELEMENTS IN IWORK. LOCIW MUST BE LARGER
C THAN 4*P.
C
C THIS VERSION OF LDA IS DATED 10/1/84.
C
C REFERENCES:
C
C V. E. KANE, R. C. WARD, AND G. J. DAVIS, "ASSESSMENT OF LINEAR
C DEPENDENCIES IN MULTIVARIATE DATA," SIAM J. SCI. STAT.
C COMPUT., TO APPEAR
C
C R. C. WARD, G. J. DAVIS, AND V. E. KANE, "AN ALGORITHM FOR
C LINEAR DEPENDENCY ANALYSIS OF MULTIVARIATE DATA,"
C ACM TOMS, SUBMITTED.
C
C THE LDA SUPPLIED CODE CONSISTS OF THE FOLLOWING SUBROUTINES:
C
C BETA, CHKEV, DEPEND, LDA, PICVAR, PRINTM, RESFLY, SCPROD.
C
C THE LDA SUPPLIED CODE CALLS THE FOLLOWING ADDITIONAL SOFTWARE:
C
C LINPACK: DPODI, DPOFA, DPOSL, DSVDC
C EISPACK: RS
C BLAS: DAXPY, DCOPY, DDOT, DNRM2, DSCAL, DROT, DROTG, DSWAP
C FORTRAN: DABS, DLOG, DMAX1, DSIGN, DSQRT, IFIX, MIN0,
C SNGL
C
C ***
INTEGER I, I1, I10, I2, I3, I4, I5, I6, I7, I8, I9, ID,
* IERR, IPRINT, J, J1, J2, J3, J4, JOB, JOBAB, KALL,
* KRES, LDL, LDX, LMAX, LOCIW, LOCW, MAX0, N, P, P2
INTEGER IWORK(LOCIW), LVAR(LDL,LMAX), IDUM(1,1)
DOUBLE PRECISION PCERR, WORK(LOCW), X(LDX,P), DUM(1,1)
IERR = 0
ID = 1
CALL PRINTM(1, DUM, ID, ID, ID, IDUM, ID, ID, ID)
WRITE (6,99998)
WRITE (6,99997) N, P, LDX, JOBAB, IPRINT
IF (LDX.GE.MAX0(N,P)) GO TO 10
WRITE (6,99996)
IERR = 1
10 IF (MIN0(N,P).GE.2) GO TO 20
WRITE (6,99995)
IERR = 1
20 IF (IPRINT.GE.1 .AND. IPRINT.LE.4) GO TO 30
WRITE (6,99994)
IERR = 1
30 JOB = JOBAB/10
KALL = JOBAB - JOB*10
IF (JOB.EQ.3) KALL = 1
IF (JOB.GE.1 .AND. JOB.LE.3) GO TO 40
WRITE (6,99993)
GO TO 150
40 IF (KALL.EQ.0 .OR. KALL.EQ.1) GO TO 50
WRITE (6,99992)
IERR = 1
50 IF (JOB.NE.1) GO TO 60
WRITE (6,99991) PCERR
IF (PCERR.GE.0.D0.AND.PCERR.LT.100.D0) GO TO 60
WRITE (6,99990)
IERR = 1
60 IF (JOB.NE.2) GO TO 70
WRITE (6,99989) P2
IF (P2.GT.0 .AND. P2.LT.P) GO TO 70
WRITE (6,99988)
IERR = 1
70 IF (JOB.NE.3) GO TO 120
WRITE (6,99987) P2, LMAX
DO 80 I=1,P2
WRITE (6,99986) (LVAR(I,J),J=1,LMAX)
80 CONTINUE
IF (P2.GT.0 .AND. P2.LT.P) GO TO 90
WRITE (6,99988)
IERR = 1
90 DO 110 I=1,P2
DO 100 J=1,LMAX
K = LVAR(I,J)
IF (K.GE.1 .AND. K.LE.P) GO TO 100
WRITE (6,99999) I, J
IERR = 1
100 CONTINUE
110 CONTINUE
GO TO 130
120 IF (KALL.EQ.1) WRITE (6,99985) KRES, LMAX
IF (KALL.NE.0) GO TO 130
IF (LMAX.EQ.1) GO TO 130
WRITE (6,99984)
IERR = 1
130 IF (KALL.EQ.1 .AND. KRES.LT.0) KRES = 2
I1 = 1 + N*P
I2 = I1 + P
I3 = I2 + P*P
I4 = I3 + P*P
I5 = I4 + P
I6 = I5 + P*P
I7 = I6 + N
I8 = I7 + P
I9 = I8 + LMAX*5
I10 = I9 + P - 1
J1 = P + 1
J2 = J1 + P
J3 = J2 + P
J4 = J3 + P - 1
WRITE (6,99983) LOCW, I10
WRITE (6,99982) LOCIW, J4
IF (I10.LE.LOCW .AND. J4.LE.LOCIW) GO TO 140
WRITE (6,99981)
IERR = 1
140 IF (IERR.NE.0) GO TO 150
WRITE (6,99980)
CALL DEPEND(X, LDX, N, P, JOB, KALL, PCERR, P2, LVAR,
* LDL, LMAX, KRES, WORK(1), WORK(I1), WORK(I2),
* WORK(I3), WORK(I4), WORK(I5), WORK(I6), WORK(I7),
* WORK(I8), WORK(I9), IWORK(1), IWORK(J1), IWORK(J2),
* IWORK(J3), IPRINT)
150 RETURN
99999 FORMAT (/6H LVAR(, I3, 1H,, I3, 22H) LESS THAN 1 OR GREAT,
* 9HER THAN P//36H **** PROGRAM ABORTING AFTER CHECKIN,
* 22HG OTHER ARGUMENTS ****)
99998 FORMAT (1X, 16HINPUT PARAMETERS/)
99997 FORMAT (4H N =, I4/4H P =, I3/6H LDX =, I4/8H JOBAB =,
* I3/9H IPRINT =, I3)
99996 FORMAT (/45H LDX NOT GREATER THAN OR EQUAL TO MAXIMUM OF ,
* 7HN AND P//38H **** PROGRAM ABORTING AFTER CHECKING ,
* 20HOTHER ARGUMENTS ****)
99995 FORMAT (/26H N OR P NOT GREATER THAN 1//14H **** PROGRAM ,
* 44HABORTING AFTER CHECKING OTHER ARGUMENTS ****)
99994 FORMAT (/27H IPRINT NOT BETWEEN 1 AND 4//13H **** PROGRAM,
* 45H ABORTING AFTER CHECKING OTHER ARGUMENTS ****)
99993 FORMAT (/34H A IN JOBAB IS NOT BETWEEN 1 AND 3//7H **** P,
* 20HROGRAM ABORTING ****)
99992 FORMAT (/25H B IN JOBAB IS NOT 0 OR 1//15H **** PROGRAM A,
* 43HBORTING AFTER CHECKING OTHER ARGUMENTS ****)
99991 FORMAT (8H PCERR =, F7.2)
99990 FORMAT (/30H PCERR NOT BETWEEN 0. AND 100.//10H **** PROG,
* 48HRAM ABORTING AFTER CHECKING OTHER ARGUMENTS ****)
99989 FORMAT (5H P2 =, I3)
99988 FORMAT (/23H P2 NOT BETWEEN 0 AND P//17H **** PROGRAM ABO,
* 41HRTING AFTER CHECKING OTHER ARGUMENTS ****)
99987 FORMAT (5H P2 =, I3/7H LMAX =, I3/16H LVAR(P2,LMAX) =)
99986 FORMAT (1X, 40I3)
99985 FORMAT (7H KRES =, I3/7H LMAX =, I3)
99984 FORMAT (/43H B IN JOBAB IS ZERO, LMAX MUST BE SET TO 1.//
* 49H **** PROGRAM ABORTING AFTER CHECKING OTHER ARGUM,
* 9HENTS ****)
99983 FORMAT (7H LOCW =, I7, 10X, 26H(ACTUAL WORK SPACE NEEDED ,
* 1H=, I7, 1H))
99982 FORMAT (8H LOCIW =, I6, 10X, 25H(ACTUAL IWORK SPACE NEEDE,
* 3HD =, I6, 1H))
99981 FORMAT (/43H USER DID NOT ALLOCATE ENOUGH WORKING SPACE//
* 27H **** PROGRAM ABORTING ****)
99980 FORMAT (////)
END
C
C --------------------------------------------------------------
C
SUBROUTINE DEPEND(X, LDX, N, P, JOB, KALL, PCERR, P2,
* LVAR, LDL, LMAX, KRES, XSAVE, W, V, CORR, XMU, R,
* WK, WK2, OUTPUT, COVD, IVAR, INT, DVAR, PVAR, IPRINT)
C ***
C SUBROUTINE DEPEND DRIVES THE COMPUTATIONS IN THE LINEAR
C DEPENDENCY ANALYSIS CALLING ON VARIOUS SUBROUTINES TO
C COMPUTE NEEDED QUANTITIES.
C ***
INTEGER I, ID, IE, IERR, IFIX, II, INFO, IOPT, IP,
* IPRINT, ISC, IT1, IT2, IT3, J, JE, JOB, JP, K, KALL,
* KE, KP1, KRES, LAST, LDL, LDX, LIST, LKNT, LMAX,
* LOC, LRES, MATZ, N, NO, P, P1, P12, P2, P2PJ
INTEGER DVAR(P), INT(P), IVAR(P), LVAR(LDL,LMAX),
* PVAR(P), IDUM(1,1)
DOUBLE PRECISION BOUND, CRES, ENU2, PCERR, PROD, PSIK,
* BIG, PSIMIN, RLKDR, RN, RNM1, RVAL, XLAMAX, DUMS,
* TEMP
DOUBLE PRECISION CORR(P,P), COVD(P), DET(2),
* OUTPUT(LMAX,5), R(P,P), V(P,P), W(P), WK(N), WK2(P),
* X(LDX,P), XMU(P), XSAVE(N,P), DUM(1,1)
DOUBLE PRECISION DLOG, DSIGN, DSQRT
REAL SNGL
C ***
C PROLOG.
C ***
ID = 1
IF (IPRINT.EQ.4) CALL PRINTM(4, X, LDX, N, P, IDUM, ID,
* ID, ID)
RN = N
RNM1 = N-1
MATZ = 0
BIG = 1.D35
C ***
C COMPUTE MEANS.
C ***
DO 20 J=1,P
XMU(J) = 0.D0
DO 10 I=1,N
XMU(J) = XMU(J) + X(I,J)
10 CONTINUE
XMU(J) = XMU(J)/RN
20 CONTINUE
IF (IPRINT.GE.2) CALL PRINTM(2, XMU, LDX, P, 1, IDUM, ID,
* ID, ID)
C ***
C COMPUTE COVARIANCE MATRIX.
C ***
DO 50 I=1,P
DO 40 J=1,P
CORR(I,J) = 0.D0
DO 30 K=1,N
CORR(I,J) = CORR(I,J) + X(K,I)*X(K,J)
30 CONTINUE
CORR(I,J) = (CORR(I,J)-RN*XMU(I)*XMU(J))/RNM1
40 CONTINUE
50 CONTINUE
IF (IPRINT.GE.3) CALL PRINTM(3, CORR, P, P, P, IDUM, ID,
* ID, ID)
C ***
C CHECK FOR CONSTANT VARIABLES
C ***
K = 0
DO 60 I=1,P
IF (CORR(I,I).NE.0.D0) GO TO 60
IDUM(1,1) = I
CALL PRINTM(36, DUM, ID, ID, ID, IDUM, 1, 1, 1)
K = 1
60 CONTINUE
IF (K.EQ.1) GO TO 600
C ***
C COMPUTE CENTERED AND SCALED X, AND SAVE.
C ***
DO 80 I=1,N
DO 70 J=1,P
X(I,J) = (X(I,J)-XMU(J))/DSQRT(CORR(J,J))
XSAVE(I,J) = X(I,J)
70 CONTINUE
80 CONTINUE
IDUM(1,1) = 1
IF (IPRINT.EQ.4) CALL PRINTM(9, X, LDX, N, P, IDUM, 1, 1,
* 1)
C ***
C SAVE COVARIANCE DIAGONALS AND COMPUTE CORRELATION MATRIX.
C ***
DO 90 I=1,P
COVD(I) = CORR(I,I)
90 CONTINUE
DO 110 I=1,P
DO 100 J=I,P
CORR(I,J) = CORR(I,J)/DSQRT(COVD(I)*COVD(J))
CORR(J,I) = CORR(I,J)
100 CONTINUE
110 CONTINUE
IF (IPRINT.GE.3) CALL PRINTM(6, CORR, P, P, P, IDUM, ID,
* ID, ID)
C ***
C COMPUTE SVD OF CENTERED, SCALED DATA MATRIX.
C ***
IOPT = 1
CALL DSVDC(X, LDX, N, P, W, R, DUM, 1, X, LDX, WK, IOPT,
* IERR)
IF (IERR.EQ.0) GO TO 120
CALL PRINTM(13, DUM, ID, ID, ID, IDUM, ID, ID, ID)
GO TO 600
120 DO 130 I=1,P
R(I,1) = W(I)
R(I,2) = (W(I)**2)/RNM1
W(I) = R(I,2)
R(I,3) = 100.D0
IF (W(I).GE..5D0) GO TO 130
R(I,3) = 100.D0*DSQRT(W(I)/(1.D0-W(I)))
130 CONTINUE
CALL PRINTM(7, R, P, P, 3, IDUM, ID, ID, ID)
TEMP = BIG
IF (R(P,1).NE.0.D0) TEMP = R(1,1)/R(P,1)
R(1,1) = TEMP
R(1,2) = BIG
IF (W(P).NE.0.D0) R(1,2) = W(1)/W(P)
CALL PRINTM(5, R, P, 1, 2, IDUM, ID, ID, ID)
IF (IPRINT.EQ.4) CALL PRINTM(8, X, LDX, P, P, IDUM, ID,
* ID, ID)
IF (JOB.NE.1) GO TO 160
C ***
C COMPUTE P2 FROM PER CENT ERROR.
C ***
ENU2 = (PCERR/100.D0)**2
BOUND = ENU2/(1.D0+ENU2)
DO 140 I=1,P
II = P + 1 - I
IF (W(II).GT.BOUND) GO TO 150
140 CONTINUE
150 P2 = I - 1
R(1,1) = PCERR
R(2,1) = BOUND
160 IF (P2.GT.0) GO TO 170
CALL PRINTM(10, DUM, ID, ID, ID, IDUM, ID, ID, ID)
GO TO 600
170 IF (P2.LT.P) GO TO 180
CALL PRINTM(11, DUM, ID, ID, ID, IDUM, ID, ID, ID)
GO TO 600
180 P1 = P - P2
CALL PRINTM(12, R, P, 2, 1, IDUM, LMAX, JOB, P2)
IF (JOB.EQ.3) GO TO 310
C ***
C COMPUTE NEARLY OPTIMAL SET OF ESTIMATED VARIABLES.
C ***
DO 210 I=1,P2
II = P + 1 - I
DO 190 J=1,P
V(I,J) = X(J,II)
190 CONTINUE
DO 200 J=1,P2
R(I,J) = 0.D0
200 CONTINUE
R(I,I) = DSQRT(RNM1*W(II))
210 CONTINUE
DO 220 I=1,P
IVAR(I) = I
220 CONTINUE
IOPT = 0
IF (IPRINT.EQ.4) CALL PRINTM(29, V, P, P2, P, IDUM, ID,
* ID, ID)
CALL PICVAR(P2, P, V, P, IOPT, IVAR, R, P, WK, WK2, RVAL)
C ***
C COMPUTE ORDERED PREDICTOR VARIABLES.
C ***
DO 240 J=1,P1
P2PJ = P2 + J
K = IVAR(P2PJ)
DO 230 I=1,P1
V(I,J) = X(K,I)
230 CONTINUE
240 CONTINUE
IOPT = 1
CALL PICVAR(P1, P1, V, P, IOPT, IVAR(P2+1), DUM, 1, WK,
* WK2, DUMS)
DUM(1,1) = RVAL
CALL PRINTM(14, DUM, 1, 1, 1, IVAR, P, P2, 1)
IF (KALL.EQ.1) GO TO 270
CALL PRINTM(34, DUM, ID, ID, ID, IDUM, ID, ID, ID)
CRES = RVAL
DO 250 I=1,P
PVAR(I) = IVAR(I)
250 CONTINUE
DO 260 I=1,P2
IT1 = P2 + 1 - I
DVAR(I) = IVAR(IT1)
260 CONTINUE
LKNT = 0
GO TO 380
C ***
C ORDER VARIABLES.
C ***
270 IF (P1.EQ.1) GO TO 290
P12 = P1/2
DO 280 I=1,P12
IT1 = P2 + I
IT2 = P + 1 - I
K = IVAR(IT1)
IVAR(IT1) = IVAR(IT2)
IVAR(IT2) = K
280 CONTINUE
C ***
C START SEQUENTIAL VARIABLE REMOVAL PROCESS.
C ***
290 DO 300 I=1,P2
INT(I) = P + 1 - I
300 CONTINUE
IDUM(1,1) = 1
CALL PRINTM(32, DUM, ID, ID, ID, IDUM, 1, 1, 1)
310 LKNT = 0
320 CONTINUE
DO 330 I=1,P2
IT1 = INT(I)
DVAR(I) = IVAR(IT1)
IF (JOB.EQ.3) DVAR(I) = LVAR(I,LKNT+1)
330 CONTINUE
DO 360 I=1,P2
II = P + 1 - I
DO 340 J=1,P
V(I,J) = X(J,II)
340 CONTINUE
DO 350 J=1,P2
R(I,J) = 0.D0
350 CONTINUE
R(I,I) = DSQRT(RNM1*W(II))
360 CONTINUE
DO 370 I=1,P
PVAR(I) = I
370 CONTINUE
IF (IPRINT.GE.3) CALL PRINTM(15, DUM, ID, ID, ID, DVAR,
* P, P2, 1)
CALL RESFLY(P2, P, V, P, PVAR, DVAR, JOB, KRES, R, P, WK,
* WK2, RVAL, LRES, IPRINT)
CRES = WK(1)
IF (LRES.EQ.0 .AND. CRES.LT.RVAL) RVAL = CRES
IF (IPRINT.EQ.4) CALL PRINTM(26, DUM, ID, ID, ID, PVAR,
* P, P, 1)
DUM(1,1) = CRES
IF (IPRINT.GE.3 .AND. LRES.EQ.0) CALL PRINTM(16, DUM, 1,
* 1, 1, IDUM, ID, ID, ID)
C ***
C CHECK EIGENVALUE INEQUALITIES FOR SMALL RESIDUALS AND PROCEED.
C ***
IF (JOB.EQ.3) GO TO 380
IF (LRES.NE.0) GO TO 540
380 DO 400 I=1,P
IE = P + 1 - I
IP = PVAR(IE)
DO 390 J=1,P
JE = P + 1 - J
JP = PVAR(JE)
V(I,J) = CORR(IP,JP)
390 CONTINUE
400 CONTINUE
IOPT = 0
CALL CHKEV(V, P, P, P1, IOPT, R, WK, WK2, PSIMIN, XLAMAX,
* PSIK, IPRINT)
C ***
C COMPUTE LIKELIHOOD RATIO.
C ***
RLKDR = 0.D0
IF (N.LE.P) GO TO 430
DO 420 I=1,P2
IE = P2 + 1 - I
IP = PVAR(IE)
DO 410 J=I,P2
JE = P2 + 1 - J
JP = PVAR(JE)
R(I,J) = V(I,J)*DSQRT(COVD(IP)*COVD(JP))
R(J,I) = R(I,J)
410 CONTINUE
WK(I) = 1.D0/R(I,I)
420 CONTINUE
IF (IPRINT.EQ.4) CALL PRINTM(35, R, P, P2, P2, IDUM, ID,
* ID, ID)
CALL SCPROD(WK, P2, PROD, ISC)
CALL DPOFA(R, P, P2, INFO)
IDUM(1,1) = INFO
IF (INFO.NE.0) CALL PRINTM(33, DUM, ID, ID, ID, IDUM, 1,
* 1, 1)
CALL DPODI(R, P, P2, DET, 10)
PROD = PROD*DET(1)
ISC = ISC + IFIX(SNGL(DET(2)+DSIGN(.1D0,DET(2))))
RLKDR = -RN*DLOG(PROD*10.D0**ISC)
430 LKNT = LKNT + 1
C ***
C UPDATE ORDERED OUTPUT ARRAYS.
C ***
LAST = LKNT - 1
IF (LKNT.GT.1) GO TO 440
LOC = 1
GO TO 510
440 IF (LKNT.LE.LMAX) GO TO 450
IF (CRES.GT.OUTPUT(LMAX,1)) GO TO 540
LAST = LMAX
450 DO 500 K=1,LAST
IF (CRES.GE.OUTPUT(K,1)) GO TO 500
KE = LKNT
IF (LKNT.GE.LMAX) KE = LMAX
IF (K.EQ.LMAX) GO TO 490
KP1 = K + 1
DO 480 I=KP1,KE
IT1 = KE + KP1 - I
IT2 = IT1 - 1
DO 460 J=1,5
OUTPUT(IT1,J) = OUTPUT(IT2,J)
460 CONTINUE
DO 470 J=1,P2
LVAR(J,IT1) = LVAR(J,IT2)
470 CONTINUE
480 CONTINUE
490 LOC = K
GO TO 510
500 CONTINUE
LOC = LKNT
510 DO 520 I=1,P2
II = P2 + 1 - I
LVAR(I,LOC) = DVAR(II)
520 CONTINUE
OUTPUT(LOC,1) = CRES
OUTPUT(LOC,2) = PSIK
TEMP = BIG
IF (XLAMAX.NE.0.D0) TEMP = PSIMIN/XLAMAX
OUTPUT(LOC,3) = TEMP
TEMP = 1.D0
IF (W(P1).NE.0.D0) TEMP = 0.D0
IF (W(P1+1).NE.0.D0) TEMP = OUTPUT(LOC,3)/(W(P1)/W(P1+1))
IF (XLAMAX.EQ.0.D0) TEMP = 1.D0
OUTPUT(LOC,4) = TEMP
OUTPUT(LOC,5) = RLKDR
IF (JOB.EQ.3) GO TO 530
IF (KALL.EQ.0) GO TO 570
GO TO 540
530 IF (LKNT.LT.LMAX) GO TO 320
GO TO 570
C ***
C DETERMINE NEXT SET OF ESTIMATED VARIABLES TO TRY.
C ***
540 DO 560 I=1,P2
IT1 = P2 + 1 - I
IF (LRES.NE.0 .AND. I.EQ.1) IT1 = LRES
IF (INT(IT1).EQ.I) GO TO 560
INT(IT1) = INT(IT1) - 1
IF (IT1.EQ.P2) GO TO 320
IT2 = IT1 + 1
DO 550 J=IT2,P2
IT3 = J + 1 - IT2
INT(J) = INT(IT1) - IT3
550 CONTINUE
GO TO 320
560 CONTINUE
IDUM(1,1) = 2
CALL PRINTM(32, DUM, ID, ID, ID, IDUM, 1, 1, 1)
570 LIST = LMAX
IF (LIST.GT.LKNT) LIST = LKNT
CALL PRINTM(17, OUTPUT, LMAX, LIST, 5, LVAR, LDL, P2,
* LIST)
C ***
C COMPUTE BETA FOR BEST LMAX RESIDUAL CASES.
C ***
DO 590 I=1,P
DO 580 J=I,P
CORR(I,J) = CORR(I,J)*DSQRT(COVD(I)*COVD(J))
CORR(J,I) = CORR(I,J)
580 CONTINUE
590 CONTINUE
IF (IPRINT.EQ.4) CALL PRINTM(3, CORR, P, P, P, IDUM, ID,
* ID, ID)
NO = LIST
CALL BETA(LVAR, LDL, NO, XSAVE, N, P, CORR, P, P2, XMU,
* X, LDX, V, P, W, WK, WK2, PVAR, IPRINT)
600 CONTINUE
RETURN
END
C
C --------------------------------------------------------------
C
SUBROUTINE PICVAR(P2, P, V, MV, IOPT, PVAR, RES, MR, WO,
* WN, RNORM)
C ***
C SUBROUTINE PICVAR COMPUTES THE BENCHMARK SOLUTION BY A QR
C FACTORIZATION. IF IOPT = 0, THE ESTIMATED VARIABLES ARE
C COMPUTED; IF IOPT = 1, THE PREDICTOR VARIABLES ARE COMPUTED.
C ***
INTEGER I, IOPT, J, JP, K, KM1, L, LEN, LP1, MAXJ, MR,
* MV, P, P2, P2P1
INTEGER PVAR(P)
DOUBLE PRECISION CNORM, FAC, GNORM, ONEPS, PROD, RNORM
DOUBLE PRECISION RES(MR,P2), V(MV,P), WN(P), WO(P)
DOUBLE PRECISION DABS, DDOT, DMAX1, DNRM2, DSIGN, DSQRT
C ***
C COMPUTE COLUMN NORMS.
C ***
DO 10 J=1,P
WN(J) = DNRM2(P2,V(1,J),1)
WO(J) = WN(J)
10 CONTINUE
C ***
C REDUCE TO TRAPEZOIDAL FORM.
C ***
DO 70 L=1,P2
C ***
C PIVOT COLUMNS.
C ***
GNORM = 0.D0
MAXJ = L
DO 20 J=L,P
IF (WN(J).LE.GNORM) GO TO 20
GNORM = WN(J)
MAXJ = J
20 CONTINUE
IF (MAXJ.EQ.L) GO TO 30
CALL DSWAP(P2, V(1,L), 1, V(1,MAXJ), 1)
WN(MAXJ) = WN(L)
WO(MAXJ) = WO(L)
JP = PVAR(MAXJ)
PVAR(MAXJ) = PVAR(L)
PVAR(L) = JP
30 CONTINUE
WN(L) = 0.D0
IF (L.EQ.P2) GO TO 70
C ***
C REDUCE COLUMN L.
C ***
LEN = P2 - L + 1
CNORM = DNRM2(LEN,V(L,L),1)
IF (CNORM.EQ.0.D0) GO TO 70
IF (V(L,L).NE.0.D0) CNORM = DSIGN(CNORM,V(L,L))
CALL DSCAL(LEN, 1.D0/CNORM, V(L,L), 1)
V(L,L) = 1.D0 + V(L,L)
C ***
C APPLY TRANSFORMATION TO OTHER COLUMNS.
C ***
LP1 = L + 1
DO 50 J=LP1,P
PROD = -DDOT(LEN,V(L,L),1,V(L,J),1)/V(L,L)
CALL DAXPY(LEN, PROD, V(L,L), 1, V(L,J), 1)
IF (WN(J).EQ.0.D0) GO TO 50
FAC = 1.D0 - (DABS(V(L,J))/WN(J))**2
FAC = DMAX1(FAC,0.D0)
ONEPS = 1.D0 + .05D0*FAC*(WN(J)/WO(J))**2
IF (ONEPS.EQ.1.D0) GO TO 40
WN(J) = WN(J)*DSQRT(FAC)
GO TO 50
40 CONTINUE
WN(J) = DNRM2(P2-L,V(L+1,J),1)
WO(J) = WN(J)
50 CONTINUE
IF (IOPT.EQ.1) GO TO 70
C ***
C APPLY TRANSFORMATION TO RESIDUAL MATRIX.
C ***
DO 60 J=1,P2
PROD = -DDOT(LEN,V(L,L),1,RES(L,J),1)/V(L,L)
CALL DAXPY(LEN, PROD, V(L,L), 1, RES(L,J), 1)
60 CONTINUE
WN(L) = V(L,L)
V(L,L) = -CNORM
70 CONTINUE
IF (IOPT.EQ.1) GO TO 140
C ***
C COMPUTE DEPENDENCIES.
C ***
P2P1 = P2 + 1
DO 90 J=P2P1,P
V(P2,J) = V(P2,J)/V(P2,P2)
DO 80 L=2,P2
K = P2 - L + 2
KM1 = K - 1
CALL DAXPY(KM1, -V(K,J), V(1,K), 1, V(1,J), 1)
V(KM1,J) = V(KM1,J)/V(KM1,KM1)
80 CONTINUE
90 CONTINUE
C ***
C APPLY INVERSE TO RESIDUAL MATRIX.
C ***
DO 110 J=1,P2
RES(P2,J) = RES(P2,J)/V(P2,P2)
DO 100 L=2,P2
K = P2 - L + 2
KM1 = K - 1
CALL DAXPY(KM1, -RES(K,J), V(1,K), 1, RES(1,J),
* 1)
RES(KM1,J) = RES(KM1,J)/V(KM1,KM1)
100 CONTINUE
110 CONTINUE
C ***
C COMPUTE RESIDUAL NORM.
C ***
RNORM = 0.D0
DO 130 I=1,P2
DO 120 J=1,P2
RNORM = RNORM + RES(I,J)**2
120 CONTINUE
130 CONTINUE
RNORM = DSQRT(RNORM)
140 CONTINUE
RETURN
END
C
C --------------------------------------------------------------
C
SUBROUTINE RESFLY(P2, P, V, MV, PVAR, DVAR, JOB, KRES,
* RES, MR, WN, VT, RMIN, LRES, IPRINT)
C ***
C SUBROUTINE RESFLY COMPUTES THE REDUCED SPACE RESIDUAL FOR A
C PARTICULAR SET OF ESTIMATED VARIABLES BY A SEQUENTIAL VARIABLE
C DELECTION PROCESS. LARGE RESIDUALS ARE DETECTED EARLY IF
C POSSIBLE.
C ***
INTEGER I, ID, IPRINT, J, JOB, JP, K, KRES, L, LEN, LL,
* LM1, LM2, LMI, LMI1, LRES, MR, MV, NJ, P, P2
INTEGER DVAR(P2), PVAR(P), IDUM(1,1)
DOUBLE PRECISION CNORM, PROD, RMIN, RNORM, TEMP, TOOBIG
DOUBLE PRECISION RES(MR,P2), V(MV,P), VT(P2), WN(P),
* DUM(1,1)
DOUBLE PRECISION DDOT, DNRM2, DSIGN, DSQRT
ID = 1
LRES = 0
IF (JOB.NE.3) TOOBIG = RMIN*(10.D0**KRES)
C ***
C REDUCE TO TRAPEZOIDAL FORM WHILE UPDATING RESIDUALS.
C ***
DO 160 L=1,P2
NJ = DVAR(L)
DO 10 LL=L,P
IF (NJ.EQ.PVAR(LL)) GO TO 20
10 CONTINUE
20 NJ = LL
IF (NJ.EQ.L) GO TO 30
C ***
C PIVOT COLUMNS.
C ***
CALL DSWAP(P2, V(1,L), 1, V(1,NJ), 1)
JP = PVAR(NJ)
PVAR(NJ) = PVAR(L)
PVAR(L) = JP
30 CONTINUE
IF (L.EQ.1) GO TO 50
C ***
C APPLY PREVIOUS TRANSFORMATIONS TO COLUMN L OF V AND RES.
C ***
LM1 = L - 1
DO 40 K=1,LM1
LEN = P2 - K + 1
TEMP = V(K,K)
V(K,K) = WN(K)
PROD = -DDOT(LEN,V(K,K),1,V(K,L),1)/V(K,K)
CALL DAXPY(LEN, PROD, V(K,K), 1, V(K,L), 1)
PROD = -DDOT(LEN,V(K,K),1,RES(K,L),1)/V(K,K)
CALL DAXPY(LEN, PROD, V(K,K), 1, RES(K,L), 1)
V(K,K) = TEMP
40 CONTINUE
50 CONTINUE
IF (L.LT.P2) GO TO 60
IF (V(L,L).EQ.0.D0) GO TO 170
GO TO 80
C ***
C REDUCE COLUMN L.
C ***
60 LEN = P2 - L + 1
CNORM = DNRM2(LEN,V(L,L),1)
IF (CNORM.EQ.0.D0) GO TO 170
IF (V(L,L).NE.0.D0) CNORM = DSIGN(CNORM,V(L,L))
CALL DSCAL(LEN, 1.D0/CNORM, V(L,L), 1)
V(L,L) = 1.D0 + V(L,L)
C ***
C APPLY TRANSFORMATION TO FIRST L COLUMNS OF RES.
C ***
DO 70 J=1,L
PROD = -DDOT(LEN,V(L,L),1,RES(L,J),1)/V(L,L)
CALL DAXPY(LEN, PROD, V(L,L), 1, RES(L,J), 1)
70 CONTINUE
WN(L) = V(L,L)
V(L,L) = -CNORM
80 CONTINUE
C ***
C COMPUTE FIRST L COMPONENTS IN COLUMN L OF RESIDUAL MATRIX.
C ***
RES(L,L) = RES(L,L)/V(L,L)
IF (L.EQ.1) GO TO 130
CALL DAXPY(LM1, -RES(L,L), V(1,L), 1, RES(1,L), 1)
RES(LM1,L) = RES(LM1,L)/V(LM1,LM1)
VT(LM1) = V(LM1,L)/V(LM1,LM1)
IF (L.EQ.2) GO TO 110
LM2 = L - 2
DO 90 I=1,LM2
VT(I) = V(I,L)
90 CONTINUE
DO 100 I=2,LM1
LMI1 = LM1 - I + 2
LMI = LMI1 - 1
CALL DAXPY(LMI, -RES(LMI1,L), V(1,LMI1), 1,
* RES(1,L), 1)
RES(LMI,L) = RES(LMI,L)/V(LMI,LMI)
CALL DAXPY(LMI, -VT(LMI1), V(1,LMI1), 1, VT(1),
* 1)
VT(LMI) = VT(LMI)/V(LMI,LMI)
100 CONTINUE
110 CONTINUE
C ***
C COMPUTE TOP LEFT L BY L-1 RESIDUAL SUBMATRIX.
C ***
CALL DSCAL(LM1, 1.D0/V(L,L), RES(L,1), MR)
DO 120 I=1,LM1
CALL DAXPY(LM1, -VT(I), RES(L,1), MR, RES(I,1),
* MR)
120 CONTINUE
130 CONTINUE
C ***
C COMPUTE SIZE OF CURRENT RESIDUAL MATRIX.
C ***
RNORM = 0.D0
DO 150 I=1,L
DO 140 J=1,L
RNORM = RNORM + RES(I,J)**2
140 CONTINUE
150 CONTINUE
RNORM = DSQRT(RNORM)
IDUM(1,1) = L
DUM(1,1) = RNORM
IF (IPRINT.GE.3) CALL PRINTM(27, DUM, 1, 1, 1, IDUM,
* 1, 1, 1)
IF (JOB.EQ.3) GO TO 160
IF (RNORM.GT.TOOBIG) GO TO 180
160 CONTINUE
WN(1) = RNORM
GO TO 190
170 IF (JOB.NE.3) GO TO 180
C ***
C THE FOLLOWING NUMBER REPRESENTS AN INFINITE RESIDUAL NORM. IT
C SHOULD NOT PRODUCE OVERFLOW ON MOST KNOWN COMPUTERS.
C
WN(1) = 1.E+30
C
C ***
CALL PRINTM(37, DUM, ID, ID, ID, DVAR, P2, L, 1)
180 LRES = L
IDUM(1,1) = L
IF (IPRINT.GE.3) CALL PRINTM(28, DUM, ID, ID, ID, IDUM,
* 1, 1, 1)
190 CONTINUE
RETURN
END
C
C --------------------------------------------------------------
C
SUBROUTINE PRINTM(KODE, A, LD, M, N, IA, LDI, IM, IN)
C ***
C SUBROUTINE PRINTM PRINTS MOST OF THE OUTPUT FROM THE LINEAR
C DEPENDENCY ANALYSIS. (THE REMAINING OUTPUT IS PRINTED IN
C SUBROUTINES LDA AND BETA.)
C ***
INTEGER I, IE, IEP1, IM, IMP, IN, ISKIP, ISYM, J, K,
* KODE, LD, LDI, M, N, NCOLE, NCOLS, P1
INTEGER IA(LDI,IN)
DOUBLE PRECISION A(LD,N), CUM, SUM, XP
DATA IBLK /1H /, IAST /1H*/
GO TO (10, 20, 40, 50, 60, 70, 80, 100, 110, 120, 130,
* 140, 170, 180, 190, 200, 210, 210, 240, 250, 260,
* 270, 280, 290, 300, 310, 320, 330, 340, 350, 360,
* 370, 380, 390, 400, 410, 420), KODE
C
10 WRITE (6,99998)
GO TO 470
C
20 WRITE (6,99997)
DO 30 I=1,M
WRITE (6,99996) A(I,1)
30 CONTINUE
WRITE (6,99995)
GO TO 470
C
40 WRITE (6,99994)
GO TO 430
C
50 WRITE (6,99993)
GO TO 430
C
60 WRITE (6,99992) A(1,1), A(1,2)
WRITE (6,99995)
GO TO 470
C
70 WRITE (6,99991)
GO TO 430
C
80 WRITE (6,99990)
XP = M
SUM = 0.D0
DO 90 I=1,M
SUM = SUM + A(I,2)
CUM = (SUM/XP)*100.D0
WRITE (6,99989) A(I,1), A(I,2), A(I,3), CUM
90 CONTINUE
WRITE (6,99995)
GO TO 470
C
100 WRITE (6,99988)
GO TO 430
C
110 IF (IA(1,1).EQ.1) WRITE (6,99987)
IF (IA(1,1).EQ.2) WRITE (6,99986)
IF (IA(1,1).EQ.3) WRITE (6,99985) IM, IN
GO TO 430
C
120 WRITE (6,99984)
GO TO 470
C
130 WRITE (6,99983)
GO TO 470
C
140 WRITE (6,99982) IN
IF (IM.NE.1) GO TO 150
WRITE (6,99981) A(1,1), A(2,1)
GO TO 470
150 IF (IM.NE.2) GO TO 160
WRITE (6,99980)
GO TO 470
160 WRITE (6,99979) LDI
GO TO 470
C
170 WRITE (6,99978)
GO TO 470
C
180 WRITE (6,99977) IM
IE = IM
IF (IM.GT.33) IE = 33
WRITE (6,99976) (IA(I,1),I=1,IE)
IF (IM.GT.33) WRITE (6,99975) (IA(I,1),I=34,IM)
IMP = IM + 1
IE = LDI
P1 = LDI - IM
IF (P1.GT.33) IE = IM + 33
WRITE (6,99974) (IA(I,1),I=IMP,IE)
IEP1 = IE + 1
IF (P1.GT.33) WRITE (6,99975) (IA(I,1),I=IEP1,LDI)
WRITE (6,99973) A(1,1)
WRITE (6,99995)
GO TO 470
C
190 IE = IM
IF (IM.GT.33) IE = 33
WRITE (6,99972) (IA(I,1),I=1,IE)
IF (IM.GT.33) WRITE (6,99975) (IA(I,1),I=34,IM)
GO TO 470
C
200 WRITE (6,99971) A(1,1)
GO TO 470
C
210 WRITE (6,99970)
ISKIP = 1
DO 230 I=1,M
ISYM = IBLK
IF (A(I,3).GT.1.D0) GO TO 220
ISYM = IAST
ISKIP = 0
220 IE = IM
IF (IM.GT.15) IE = 15
WRITE (6,99969) (A(I,J),J=1,3), ISYM,
* (A(I,J),J=4,5), (IA(J,I),J=1,IE)
IF (IM.GT.15) WRITE (6,99967) (IA(J,I),J=16,IM)
230 CONTINUE
IF (ISKIP.EQ.0) WRITE (6,99968)
WRITE (6,99966)
GO TO 470
C
240 WRITE (6,99965) IA(1,1)
GO TO 470
C
250 WRITE (6,99964) IA(1,1)
GO TO 470
C
260 WRITE (6,99963)
GO TO 430
C
270 WRITE (6,99962)
GO TO 430
C
280 WRITE (6,99961)
GO TO 430
C
290 WRITE (6,99960)
WRITE (6,99996) (A(I,1),I=1,M)
GO TO 470
C
300 WRITE (6,99959)
WRITE (6,99996) (A(I,1),I=1,M)
GO TO 470
C
310 IE = IM
IF (IM.GT.30) IE = 30
WRITE (6,99958) (IA(I,1),I=1,IE)
IF (IM.GT.30) WRITE (6,99957) (IA(I,1),I=31,IM)
GO TO 470
C
320 WRITE (6,99956) IA(1,1), A(1,1)
GO TO 470
C
330 WRITE (6,99955) IA(1,1)
GO TO 470
C
340 WRITE (6,99954)
GO TO 430
C
350 IF (IA(1,1).EQ.1) WRITE (6,99953)
IF (IA(1,1).EQ.2) WRITE (6,99952)
GO TO 430
C
360 WRITE (6,99951) A(1,1)
GO TO 470
C
370 IF (IA(1,1).EQ.1) WRITE (6,99950)
IF (IA(1,1).EQ.2) WRITE (6,99949)
GO TO 470
C
380 WRITE (6,99948) IA(1,1)
GO TO 470
C
390 WRITE (6,99947)
GO TO 470
C
400 WRITE (6,99946)
GO TO 430
C
410 WRITE (6,99999) IA(1,1)
GO TO 470
C
420 WRITE (6,99943)
WRITE (6,99975) (IA(I,1),I=1,IM)
GO TO 470
C ***
C MATRIX PRINT UTILITY.
C ***
430 DO 450 K=1,N,4
NCOLS = K
NCOLE = K + 3
IF (NCOLE.GT.N) NCOLE = N
WRITE (6,99945) NCOLS, NCOLE
DO 440 I=1,M
WRITE (6,99944) (A(I,J),J=NCOLS,NCOLE)
440 CONTINUE
IF (NCOLE.EQ.N) GO TO 460
450 CONTINUE
460 WRITE (6,99995)
470 RETURN
99999 FORMAT (//9H VARIABLE, I3, 27H IS CONSTANT FOR ALL OBSER,
* 7HVATIONS, 38H *** CODE TERMINATES AFTER CHECKING OT,
* 13HHER VARIABLES)
99998 FORMAT (///1H1, 20X, 34H*** LINEAR DEPENDENCY ANALYSIS ***
* ///)
99997 FORMAT (28H THE VECTOR OF MEANS IS...../)
99996 FORMAT (1PD16.6)
99995 FORMAT (//)
99994 FORMAT (30H THE COVARIANCE MATRIX IS.....)
99993 FORMAT (33H THE ORIGINAL DATA MATRIX IS.....)
99992 FORMAT (19H CONDITION NUMBERS., 4X, 18HCENTERED, SCALED X,
* 12H MATRIX....., 1PD12.5/30X, 19HCORRELATION MATRIX.,
* 4H...., 1PD12.5/)
99991 FORMAT (31H THE CORRELATION MATRIX IS.....)
99990 FORMAT (1X, 18HSINGULAR VALUES OF, 8X, 14HEIGENVALUES OF,
* 7X, 14HPERCENT ERRORS, 4X, 21HEIGENVALUE CUMULATIVE/
* 1X, 18HCENTERED, SCALED X, 6X, 18HCORRELATION MATRIX,
* 4X, 16H(DIAGONAL MODEL), 9X, 10HPERCENTAGE/)
99989 FORMAT (2X, 1PD14.6, 10X, 1PD14.6, 12X, 0PF6.1, 15X, F6.1)
99988 FORMAT (46H SINGULAR VECTORS OF CORRELATION MATRIX ARE...,
* 2H..)
99987 FORMAT (/44H THE CENTERED AND SCALED DATA MATRIX IS.....)
99986 FORMAT (/45H THE PERMUTED CENTERED AND SCALED DATA MATRIX,
* 8H IS.....)
99985 FORMAT (/26H THE DATA ARRAY CONTAINING/4X, 10H THE FIRST,
* I4, 32H COLUMNS OF PREDICTOR VARIABLES,/4X, 6H THE L,
* 4HAST , I4, 22H COLUMNS OF RESIDUALS.)
99984 FORMAT (46H ERROR. COMPUTED P2 LESS THAN OR EQUAL TO ZER,
* 2HO.)
99983 FORMAT (46H ERROR. COMPUTED P2 GREATER THAN OR EQUAL TO ,
* 2HP.)
99982 FORMAT (/30H NUMBER OF DEPENDENCIES (P2) =, I3)
99981 FORMAT (6X, 29HDETERMINED FROM USER SUPPLIED, 9H ALLOWABL,
* 21HE PERCENTAGE ERROR OF, F6.1/6X, 8HYIELDING,
* 34H A BOUND ON THE SINGULAR VALUES OF, 1PD14.6/)
99980 FORMAT (6X, 17HSPECIFIED BY USER/)
99979 FORMAT (6X, 28HSPECIFIED BY USER ALONG WITH, I3, 6H SETS ,
* 30HOF ESTIMATED VARIABLES TO TEST/)
99978 FORMAT (18H ERROR FROM DSVDC./)
99977 FORMAT (/28H BENCHMARK SOLUTION FOR P2 =, I3//6X,
* 13HTHE FOLLOWING, 31H TWO SETS OF VARIABLES ARE ORDE,
* 4HRED./6X, 15HTHE MOST LIKELY, 18H ONES TO BE IN THA,
* 23HT SET ARE LISTED FIRST./)
99976 FORMAT (6X, 23HP2 ESTIMATED VARIABLES-, 2X, 33I3)
99975 FORMAT (31X, 33I3)
99974 FORMAT (6X, 23HP1 PREDICTOR VARIABLES-, 2X, 33I3)
99973 FORMAT (/6X, 41HAPPROXIMATE (REDUCED SPACE) RESIDUAL IS..,
* 1H., 2H.., 1PD16.6)
99972 FORMAT (//29H***** FOR ESTIMATED VARIABLES, 2X, 33I3/)
99971 FORMAT (31H REDUCED SPACE RESIDUAL IS....., 1PD16.8/)
99970 FORMAT (/////1X, 13HREDUCED SPACE, 4X, 13HPSI CONDITION,
* 22X, 10HPSI/LAMBDA/3X, 8HRESIDUAL, 11X, 6HNUMBER,
* 6X, 14HPSI/LAMBDA GAP, 3X, 14HNORMALIZED GAP, 2X,
* 16HLIKELIHOOD RATIO, 3X, 19HESTIMATED VARIABLES/)
99969 FORMAT (1X, 2(1PD12.5, 5X), 1PD12.5, 1X, A1, 3X,
* 2(1PD12.5, 5X), 15I3)
99968 FORMAT (/45H * DENOTES CASES WHERE LARGEST EIGENVALUE OF,
* 1H , 42HLAMBDA EXCEEDS SMALLEST EIGENVALUE OF PSI.)
99967 FORMAT (86X, 15I3)
99966 FORMAT (///)
99965 FORMAT (46H MATRIX NOT POSITIVE DEFINITE IN CHKEV. IERR ,
* 1H=, I4)
99964 FORMAT (44H ERROR RETURN FROM EISPACK IN CHKEV. IERR =,
* I4)
99963 FORMAT (23H THE PSI MATRIX IS.....)
99962 FORMAT (26H THE LAMBDA MATRIX IS.....)
99961 FORMAT (46H THE PARTIALLY FACTORIZED SIGMA MATRIX IS.....)
99960 FORMAT (35H THE EIGENVALUES OF LAMBDA ARE...../)
99959 FORMAT (32H THE EIGENVALUES OF PSI ARE...../)
99958 FORMAT (37H VARIABLE ORDERING OUT OF RESFLY....., 30I3)
99957 FORMAT (37X, 30I3)
99956 FORMAT (7H STEP =, I3, 10X, 10HRESIDUAL =, 1PD15.7)
99955 FORMAT (32H LARGE RESIDUAL DETECTED AT STEP, I3)
99954 FORMAT (46H THE SINGULAR VECTORS ENTERING PICVAR ARE.....)
99953 FORMAT (40H THE PERMUTED CORRELATION MATRIX IS.....)
99952 FORMAT (39H THE PERMUTED COVARIANCE MATRIX IS.....)
99951 FORMAT (//17H RESIDUAL IS....., 1PD16.8/)
99950 FORMAT (///43H BEGIN EXAMINING ALL POSSIBLE SUBSETS OF P2,
* 10H VARIABLES/)
99949 FORMAT (/46H ALL POSSIBLE SUBSETS OF P2 VARIABLES EXAMINED
* ///)
99948 FORMAT (/8H ERROR =, I3, 3X, 25HFROM DPOFA IN LIKELIHOOD ,
* 6HRATIO , 11HCOMPUTATION)
99947 FORMAT (/45H ONLY BENCHMARK SOLUTION REQUESTED. ALL POSSI,
* 3HBLE, 21H SUBSETS NOT EXAMINED///)
99946 FORMAT (/35H THE UNSCALED LAMBDA MATRIX IS.....)
99945 FORMAT (/8H COLUMNS, I5, 8H THROUGH, I5, 5H...../)
99944 FORMAT (4(1PD16.6))
99943 FORMAT (46H SINGULAR V22 RESULTS FROM ESTIMATED VARIABLES,
* 5H.....)
END
C
C --------------------------------------------------------------
C
SUBROUTINE CHKEV(SIG, LDS, P, P1, IOPT, W, TEM1, TEM2,
* PSIMIN, XLAMAX, PSIK, IPRINT)
C ***
C SUBROUTINE CHKEV COMPUTES THE PSI AND LAMBDA MATRICES FOR A
C PARTICULAR SET OF ESTIMATED VARIABLES. THE MINIMUM EIGENVALUE
C OF PSI AND THE MAXIMUM EIGENVALUE OF LAMBDA ARE COMPUTED.
C ***
INTEGER I, ID, IERR, IOPT, IP1, IPRINT, J, JM1, JP1, K,
* LDS, LEN, MATZ, P, P1, P2
INTEGER IDUM(1,1)
DOUBLE PRECISION PSIK, PSIMIN, S, T, XLAMAX
DOUBLE PRECISION SIG(LDS,P), TEM1(P), TEM2(P), W(P),
* DUM(1,1)
DOUBLE PRECISION DDOT, DSQRT
ID = 1
P2 = P - P1
IDUM(1,1) = 1
IF (IPRINT.EQ.4 .AND. IOPT.EQ.0) CALL PRINTM(30, SIG,
* LDS, P, P, IDUM, 1, 1, 1)
IDUM(1,1) = 2
IF (IPRINT.EQ.4 .AND. IOPT.EQ.1) CALL PRINTM(30, SIG,
* LDS, P, P, IDUM, 1, 1, 1)
MATZ = 0
C ***
C COMPUTE PSI (CORRELATION MATRIX OF ORDER P1) AND
C LAMBDA (SCHUR COMPLEMENT OF PSI).
C ***
DO 30 J=1,P
IERR = J
S = 0.D0
JM1 = J - 1
IF (JM1.LT.1) GO TO 20
DO 10 K=1,JM1
LEN = K - 1
IF (K.GT.P1) LEN = P1
T = SIG(K,J) - DDOT(LEN,SIG(1,K),1,SIG(1,J),1)
IF (K.LE.P1) T = T/SIG(K,K)
SIG(K,J) = T
IF (K.LE.P1) S = S + T*T
10 CONTINUE
20 S = SIG(J,J) - S
SIG(J,J) = S
IF (J.GT.P1) GO TO 30
IF (S.LE.0.D0) GO TO 40
SIG(J,J) = DSQRT(S)
30 CONTINUE
IERR = 0
IF (IPRINT.EQ.4) CALL PRINTM(23, SIG, LDS, P, P, IDUM,
* ID, ID, ID)
40 IF (IERR.EQ.0) GO TO 50
IDUM(1,1) = IERR
CALL PRINTM(19, DUM, ID, ID, ID, IDUM, 1, 1, 1)
GO TO 110
50 IF (IOPT.EQ.1) GO TO 110
C ***
C COMPUTE EIGENVALUESS OF PSI.
C ***
DO 60 I=1,P1
SIG(I,I) = 1.D0
60 CONTINUE
IF (IPRINT.EQ.4) CALL PRINTM(21, SIG, LDS, P1, P1, IDUM,
* ID, ID, ID)
CALL RS(LDS, P1, SIG, W, MATZ, DUM, TEM1, TEM2, IERR)
IF (IERR.EQ.0) GO TO 70
IDUM(1,1) = IERR
CALL PRINTM(20, DUM, ID, ID, ID, IDUM, 1, 1, 1)
GO TO 110
70 PSIMIN = W(1)
PSIK = W(P1)/W(1)
IF (IPRINT.EQ.4) CALL PRINTM(25, W, P1, P1, 1, IDUM, ID,
* ID, ID)
C ***
C COMPUTE EIGENVALUES OF LAMBDA.
C ***
DO 90 I=1,P2
IP1 = I + P1
DO 80 J=I,P2
JP1 = J + P1
SIG(J,I) = SIG(IP1,JP1)
SIG(I,J) = SIG(J,I)
80 CONTINUE
90 CONTINUE
IF (IPRINT.EQ.4) CALL PRINTM(22, SIG, LDS, P2, P2, IDUM,
* ID, ID, ID)
CALL RS(LDS, P2, SIG, W, MATZ, DUM, TEM1, TEM2, IERR)
IF (IERR.EQ.0) GO TO 100
IDUM(1,1) = IERR
CALL PRINTM(20, DUM, ID, ID, ID, IDUM, 1, 1, 1)
GO TO 110
100 XLAMAX = W(P2)
IF (IPRINT.EQ.4) CALL PRINTM(24, W, P2, P2, 1, IDUM, ID,
* ID, ID)
110 RETURN
END
C
C --------------------------------------------------------------
C
SUBROUTINE BETA(LVAR, LDL, NO, XSAVE, N, P, CORR, LDC,
* P2, XMU, X, LDX, V, LDV, W, WK, WK2, PVAR, IPRINT)
C ***
C SUBROUTINE BETA COMPUTES AND OUTPUTS FOR THE SELECTED SETS OF
C ESTIMATED VARIABLES THE BETA MATRIX OF COEFFICIENTS, THE BETA
C VECTOR OF INTERCEPTS, 95 PER CENT CONFIDENCE INTERVALS, THE
C ZERO/NONZERO STRUCTURE OF THE BETAS, AND THE FULL SPACE
C RESIDUAL.
C ***
INTEGER I, ID, IE, IERR, IL, IOPT, IP, IPRINT, IS, J,
* JIS, JJ, JL, JN, JP, JP1, K, KK, KP1, KT, L, LDC,
* LDL, LDV, LDX, LEN, N, NO, P, P1, P2
INTEGER LSTR(36), LVAR(LDL,NO), PVAR(P), IDUM(1,1)
DOUBLE PRECISION CVAR, QUAD, RESID, RN, SRTN, SUM, DUM1,
* DUM2, DUM3
DOUBLE PRECISION CORR(LDC,P), V(LDV,P), W(P), WK(N),
* WK2(P), X(LDX,P), XMU(P), XSAVE(N,P), DUM(1,1)
DOUBLE PRECISION DABS, DSQRT
DATA IZERO /1H0/, IAST /1H*/
ID = 1
P1 = P - P2
C ***
DO 350 K=1,NO
WRITE (6,99999)
C ***
C DETERMINE VARIABLE ORDER.
C ***
JP = 0
DO 20 J=1,P
DO 10 JJ=1,P2
IF (J.EQ.LVAR(JJ,K)) GO TO 20
10 CONTINUE
JP = JP + 1
PVAR(JP) = J
20 CONTINUE
DO 30 J=1,P2
JL = P1 + J
JP = LVAR(J,K)
PVAR(JL) = JP
30 CONTINUE
C ***
C COPY PERMUTED COVARIANCE MATRIX.
C ***
DO 50 I=1,P
IP = PVAR(I)
DO 40 J=1,P
JP = PVAR(J)
V(I,J) = CORR(IP,JP)
40 CONTINUE
50 CONTINUE
C ***
C COMPUTE LAMBDA MATRIX DIAGONALS.
C ***
IOPT = 1
CALL CHKEV(V, LDV, P, P1, IOPT, W, WK2, WK, DUM1,
* DUM2, DUM3, IPRINT)
DO 60 I=1,P2
IP = P1 + I
WK(I) = V(IP,IP)
60 CONTINUE
C ***
C COMPUTE BETA BY PSI(INV)*(PSI*BETA).
C ***
DO 80 I=1,P
IP = PVAR(I)
DO 70 J=1,P
JP = PVAR(J)
V(I,J) = CORR(IP,JP)
70 CONTINUE
80 CONTINUE
CALL DPOFA(V, LDV, P1, IERR)
IF (IERR.NE.0) WRITE (6,99998) IERR
DO 90 J=1,P2
JP = P1 + J
CALL DPOSL(V, LDV, P1, V(1,JP))
90 CONTINUE
C ***
C COMPUTE QUADRATIC FORM FOR CONFIDENCE INTERVALS.
C ***
DO 100 I=1,P1
IP = PVAR(I)
W(I) = XMU(IP)
100 CONTINUE
CALL DPOSL(V, LDV, P1, W)
QUAD = 0.D0
DO 110 I=1,P1
IP = PVAR(I)
QUAD = QUAD + XMU(IP)*W(I)
110 CONTINUE
QUAD = QUAD + 1.D0
C ***
C COMPUTE DIAGONALS OF THE PSI INVERSE MATRIX.
C ***
DO 130 I=1,P1
DO 120 J=1,P1
W(J) = 0.D0
120 CONTINUE
W(I) = 1.D0
CALL DPOSL(V, LDV, P1, W)
WK2(I) = W(I)
130 CONTINUE
C ***
C COMPUTE CONFIDENCE INTERVALS FOR BETA.
C ***
RN = N
SRTN = DSQRT(RN)
DO 150 I=1,P2
IP = I + P1
DO 140 J=1,P1
V(IP,J) = 1.96D0*DSQRT(WK2(J)*WK(I))/SRTN
140 CONTINUE
150 CONTINUE
C ***
C COMPUTE CONFIDENCE INTERVALS FOR BETA0.
C ***
DO 160 I=1,P2
IP = I + P1
W(I) = 1.96D0*DSQRT(QUAD*WK(I))/SRTN
160 CONTINUE
C ***
C COMPUTE BETA0.
C ***
DO 180 I=1,P2
SUM = 0.D0
IL = I + P1
DO 170 J=1,P1
JP = PVAR(J)
SUM = SUM + V(J,IL)*XMU(JP)
170 CONTINUE
IP = PVAR(IL)
WK2(I) = XMU(IP) - SUM
180 CONTINUE
C ***
C OUTPUT RESULTS.
C ***
WRITE (6,99997)
KT = (P1+7)/8
DO 220 KK=1,KT
IS = (KK-1)*8 + 1
IE = KK*8
IF (IE.GT.P1) IE = P1
IF (KK.GT.1) GO TO 200
WRITE (6,99996) (PVAR(I),I=IS,IE)
DO 190 I=1,P2
J = I + P1
WRITE (6,99995) PVAR(J), WK2(I),
* (V(L,J),L=IS,IE)
WRITE (6,99994) W(I), (V(J,L),L=IS,IE)
190 CONTINUE
GO TO 220
200 WRITE (6,99993) (PVAR(I),I=IS,IE)
DO 210 I=1,P2
J = I + P1
WRITE (6,99992) PVAR(J), (V(L,J),L=IS,IE)
WRITE (6,99991) (V(J,L),L=IS,IE)
210 CONTINUE
220 CONTINUE
WRITE (6,99990)
KT = (P1+34)/35
DO 280 KK=1,KT
IS = (KK-1)*35 + 1
IE = KK*35
IF (IE.GT.P1) IE = P1
IF (KK.GT.1) GO TO 250
WRITE (6,99989) (PVAR(I),I=IS,IE)
LEN = IE
KP1 = LEN + 1
DO 240 I=1,P2
IP = I + P1
LSTR(1) = IZERO
IF (DABS(WK2(I)).GT.W(I)) LSTR(1) = IAST
DO 230 J=1,LEN
JP1 = J + 1
LSTR(JP1) = IZERO
IF (DABS(V(J,IP)).GT.V(IP,J))
* LSTR(JP1) = IAST
230 CONTINUE
WRITE (6,99988) PVAR(IP), (LSTR(J),J=1,KP1)
240 CONTINUE
GO TO 280
250 WRITE (6,99987) (PVAR(I),I=IS,IE)
LEN = IE - IS + 1
DO 270 I=1,P2
IP = I + P1
DO 260 J=1,LEN
JIS = J + IS - 1
LSTR(J) = IZERO
IF (DABS(V(JIS,IP)).GT.V(IP,JIS))
* LSTR(J) = IAST
260 CONTINUE
WRITE (6,99986) PVAR(IP), (LSTR(J),J=1,LEN)
270 CONTINUE
280 CONTINUE
C ***
C COMPUTE CENTERED AND SCALED BETA MATRIX.
C ***
DO 300 I=1,P1
IP = PVAR(I)
CVAR = DSQRT(CORR(IP,IP))
DO 290 J=1,P2
JN = P1 + J
JP = PVAR(JN)
V(I,JN) = V(I,JN)*(CVAR/DSQRT(CORR(JP,JP)))
290 CONTINUE
300 CONTINUE
C ***
C PERMUTE COLUMNS OF CENTERED AND SCALED X.
C ***
DO 310 J=1,P
JP = PVAR(J)
CALL DCOPY(N, XSAVE(1,JP), 1, X(1,J), 1)
310 CONTINUE
IDUM(1,1) = 2
IF (IPRINT.EQ.4) CALL PRINTM(9, X, LDX, N, P, IDUM,
* 1, 1, 1)
C ***
C COMPUTE RESIDUAL (FULL SPACE RESIDUAL).
C ***
RESID = 0.D0
DO 340 I=1,N
DO 330 J=1,P2
JP = J + P1
SUM = 0.D0
DO 320 L=1,P1
SUM = SUM + X(I,L)*V(L,JP)
320 CONTINUE
X(I,JP) = X(I,JP) - SUM
RESID = RESID + X(I,JP)**2
330 CONTINUE
340 CONTINUE
IDUM(1,1) = 3
IF (IPRINT.EQ.4) CALL PRINTM(9, X, LDX, N, P, IDUM,
* 1, P1, P2)
RESID = DSQRT(RESID)
DUM(1,1) = RESID
CALL PRINTM(31, DUM, 1, 1, 1, IDUM, ID, ID, ID)
350 CONTINUE
RETURN
99999 FORMAT (///1X, 120(1H-))
99998 FORMAT (27H IERR FROM DPOFA IN BETA IS, I3)
99997 FORMAT (///43H BETA0 - BETA MATRIX / 95 PER CENT CONFIDEN,
* 3HCE , 9HINTERVALS//)
99996 FORMAT (25X, 9HPREDICTOR/1X, 9HESTIMATED, 15X, 8HVARIABLE,
* 5HS..../1X, 9HVARIABLES, 3X, 9HINTERCEPT, 1X, 8(6X,
* I3, 4X))
99995 FORMAT (/4X, I3, 3X, 1P9D13.4)
99994 FORMAT (10X, 1P9D13.4)
99993 FORMAT (//25X, 9HPREDICTOR/1X, 9HESTIMATED, 15X, 6HVARIAB,
* 3HLES, 2X, 14H(CONTINUED).../1X, 9HVARIABLES, 13X,
* 8(6X, I3, 4X))
99992 FORMAT (/4X, I3, 16X, 1P8D13.4)
99991 FORMAT (23X, 1P8D13.4)
99990 FORMAT (///1X, 39HDEPENDENCY STRUCTURE AT 95 PER CENT CON,
* 8HFIDENCE , 5HLEVEL/11X, 24H(* = NONZERO, 0 = ZERO)/
* )
99989 FORMAT (13X, 9HPREDICTOR/1X, 9HESTIMATED, 3X, 9HVARIABLES,
* 4H..../1X, 9HVARIABLES, 3H I, 35(I3))
99988 FORMAT (/4X, I3, 3X, 36(2X, A1))
99987 FORMAT (13X, 9HPREDICTOR/1X, 9HESTIMATED, 3X, 9HVARIABLES,
* 1H , 14H(CONTINUED).../1X, 9HVARIABLES, 3X, 35I3)
99986 FORMAT (/10X, 35(2X, A1))
END
C
C --------------------------------------------------------------
C
SUBROUTINE SCPROD(X, N, PROD, ISC)
C ***
C SUBROUTINE SCPROD COMPUTES THE PRODUCT OF A SEQUENCE OF NUMBERS
C USING A SCALING PROCEDURE WHICH COMPUTES THE ANSWER AS A
C MANTISSA AND EXPONENT.
C ***
INTEGER I, ISC, N
DOUBLE PRECISION PROD, SC, X(N)
DOUBLE PRECISION DABS
PROD = 1.D0
ISC = 0
SC = 10.D0
DO 30 I=1,N
PROD = X(I)*PROD
IF (PROD.EQ.0.D0) GO TO 40
10 IF (DABS(PROD).GE.1.D0) GO TO 20
PROD = SC*PROD
ISC = ISC - 1
GO TO 10
20 IF (DABS(PROD).LT.SC) GO TO 30
PROD = PROD/SC
ISC = ISC + 1
GO TO 20
30 CONTINUE
40 CONTINUE
RETURN
END
SUBROUTINE DPODI(A,LDA,N,DET,JOB)
00000010
INTEGER LDA,N,JOB
DOUBLE PRECISION A(LDA,1)
DOUBLE PRECISION DET(2)
C
C DPODI COMPUTES THE DETERMINANT AND INVERSE OF A CERTAIN
C DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE MATRIX (SEE BELOW)
C USING THE FACTORS COMPUTED BY DPOCO, DPOFA OR DQRDC.
C
C ON ENTRY
C
C A DOUBLE PRECISION(LDA, N)
C THE OUTPUT A FROM DPOCO OR DPOFA
C OR THE OUTPUT X FROM DQRDC.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C JOB INTEGER
C = 11 BOTH DETERMINANT AND INVERSE.
C = 01 INVERSE ONLY.
C = 10 DETERMINANT ONLY.
C
C ON RETURN
C
C A IF DPOCO OR DPOFA WAS USED TO FACTOR A THEN
C DPODI PRODUCES THE UPPER HALF OF INVERSE(A) .
C IF DQRDC WAS USED TO DECOMPOSE X THEN
C DPODI PRODUCES THE UPPER HALF OF INVERSE(TRANS(X)*X)
C WHERE TRANS(X) IS THE TRANSPOSE.
C ELEMENTS OF A BELOW THE DIAGONAL ARE UNCHANGED.
C IF THE UNITS DIGIT OF JOB IS ZERO, A IS UNCHANGED.
C
C DET DOUBLE PRECISION(2)
C DETERMINANT OF A OR OF TRANS(X)*X IF REQUESTED.
C OTHERWISE NOT REFERENCED.
C DETERMINANT = DET(1) * 10.0**DET(2)
C WITH 1.0 .LE. DET(1) .LT. 10.0
C OR DET(1) .EQ. 0.0 .
C
C ERROR CONDITION
C
C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED.
C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY
C AND IF DPOCO OR DPOFA HAS SET INFO .EQ. 0 .
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C BLAS DAXPY,DSCAL
C FORTRAN MOD
C
C INTERNAL VARIABLES
C
DOUBLE PRECISION T
DOUBLE PRECISION S
INTEGER I,J,JM1,K,KP1
C
C COMPUTE DETERMINANT
C
IF (JOB/10 .EQ. 0) GO TO 70
DET(1) = 1.0D0
DET(2) = 0.0D0
S = 10.0D0
DO 50 I = 1, N
DET(1) = A(I,I)**2*DET(1)
C ...EXIT
IF (DET(1) .EQ. 0.0D0) GO TO 60
10 IF (DET(1) .GE. 1.0D0) GO TO 20
DET(1) = S*DET(1)
DET(2) = DET(2) - 1.0D0
GO TO 10
20 CONTINUE
30 IF (DET(1) .LT. S) GO TO 40
DET(1) = DET(1)/S
DET(2) = DET(2) + 1.0D0
GO TO 30
40 CONTINUE
50 CONTINUE
60 CONTINUE
70 CONTINUE
C
C COMPUTE INVERSE(R)
C
IF (MOD(JOB,10) .EQ. 0) GO TO 140
DO 100 K = 1, N
A(K,K) = 1.0D0/A(K,K)
T = -A(K,K)
CALL DSCAL(K-1,T,A(1,K),1)
KP1 = K + 1
IF (N .LT. KP1) GO TO 90
DO 80 J = KP1, N
T = A(K,J)
A(K,J) = 0.0D0
CALL DAXPY(K,T,A(1,K),1,A(1,J),1)
80 CONTINUE
90 CONTINUE
100 CONTINUE
C
C FORM INVERSE(R) * TRANS(INVERSE(R))
C
DO 130 J = 1, N
JM1 = J - 1
IF (JM1 .LT. 1) GO TO 120
DO 110 K = 1, JM1
T = A(K,J)
CALL DAXPY(K,T,A(1,J),1,A(1,K),1)
110 CONTINUE
120 CONTINUE
T = A(J,J)
CALL DSCAL(J,T,A(1,J),1)
130 CONTINUE
140 CONTINUE
RETURN
END
SUBROUTINE DPOFA(A,LDA,N,INFO)
00000010
INTEGER LDA,N,INFO
DOUBLE PRECISION A(LDA,1)
C
C DPOFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE
C MATRIX.
C
C DPOFA IS USUALLY CALLED BY DPOCO, BUT IT CAN BE CALLED
C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
C (TIME FOR DPOCO) = (1 + 18/N)*(TIME FOR DPOFA) .
C
C ON ENTRY
C
C A DOUBLE PRECISION(LDA, N)
C THE SYMMETRIC MATRIX TO BE FACTORED. ONLY THE
C DIAGONAL AND UPPER TRIANGLE ARE USED.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C ON RETURN
C
C A AN UPPER TRIANGULAR MATRIX R SO THAT A = TRANS(R)*R
C WHERE TRANS(R) IS THE TRANSPOSE.
C THE STRICT LOWER TRIANGLE IS UNALTERED.
C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
C
C INFO INTEGER
C = 0 FOR NORMAL RETURN.
C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR
C OF ORDER K IS NOT POSITIVE DEFINITE.
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C BLAS DDOT
C FORTRAN DSQRT
C
C INTERNAL VARIABLES
C
DOUBLE PRECISION DDOT,T
DOUBLE PRECISION S
INTEGER J,JM1,K
C BEGIN BLOCK WITH ...EXITS TO 40
C
C
DO 30 J = 1, N
INFO = J
S = 0.0D0
JM1 = J - 1
IF (JM1 .LT. 1) GO TO 20
DO 10 K = 1, JM1
T = A(K,J) - DDOT(K-1,A(1,K),1,A(1,J),1)
T = T/A(K,K)
A(K,J) = T
S = S + T*T
10 CONTINUE
20 CONTINUE
S = A(J,J) - S
C ......EXIT
IF (S .LE. 0.0D0) GO TO 40
A(J,J) = DSQRT(S)
30 CONTINUE
INFO = 0
40 CONTINUE
RETURN
END
SUBROUTINE DPOSL(A,LDA,N,B)
00000010
INTEGER LDA,N
DOUBLE PRECISION A(LDA,1),B(1)
C
C DPOSL SOLVES THE DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE
C SYSTEM A * X = B
C USING THE FACTORS COMPUTED BY DPOCO OR DPOFA.
C
C ON ENTRY
C
C A DOUBLE PRECISION(LDA, N)
C THE OUTPUT FROM DPOCO OR DPOFA.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C B DOUBLE PRECISION(N)
C THE RIGHT HAND SIDE VECTOR.
C
C ON RETURN
C
C B THE SOLUTION VECTOR X .
C
C ERROR CONDITION
C
C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
C A ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES
C SINGULARITY BUT IT IS USUALLY CAUSED BY IMPROPER SUBROUTINE
C ARGUMENTS. IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED
C CORRECTLY AND INFO .EQ. 0 .
C
C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
C WITH P COLUMNS
C CALL DPOCO(A,LDA,N,RCOND,Z,INFO)
C IF (RCOND IS TOO SMALL .OR. INFO .NE. 0) GO TO ...
C DO 10 J = 1, P
C CALL DPOSL(A,LDA,N,C(1,J))
C 10 CONTINUE
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C BLAS DAXPY,DDOT
C
C INTERNAL VARIABLES
C
DOUBLE PRECISION DDOT,T
INTEGER K,KB
C
C SOLVE TRANS(R)*Y = B
C
DO 10 K = 1, N
T = DDOT(K-1,A(1,K),1,B(1),1)
B(K) = (B(K) - T)/A(K,K)
10 CONTINUE
C
C SOLVE R*X = Y
C
DO 20 KB = 1, N
K = N + 1 - KB
B(K) = B(K)/A(K,K)
T = -B(K)
CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
20 CONTINUE
RETURN
END
SUBROUTINE DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
00000010
INTEGER LDX,N,P,LDU,LDV,JOB,INFO
DOUBLE PRECISION X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1)
C
C
C DSVDC IS A SUBROUTINE TO REDUCE A DOUBLE PRECISION NXP MATRIX X
C BY ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM. THE
C DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X. THE
C COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
C AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
C
C ON ENTRY
C
C X DOUBLE PRECISION(LDX,P), WHERE LDX.GE.N.
C X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
C DECOMPOSITION IS TO BE COMPUTED. X IS
C DESTROYED BY DSVDC.
C
C LDX INTEGER.
C LDX IS THE LEADING DIMENSION OF THE ARRAY X.
C
C N INTEGER.
C N IS THE NUMBER OF ROWS OF THE MATRIX X.
C
C P INTEGER.
C P IS THE NUMBER OF COLUMNS OF THE MATRIX X.
C
C LDU INTEGER.
C LDU IS THE LEADING DIMENSION OF THE ARRAY U.
C (SEE BELOW).
C
C LDV INTEGER.
C LDV IS THE LEADING DIMENSION OF THE ARRAY V.
C (SEE BELOW).
C
C WORK DOUBLE PRECISION(N).
C WORK IS A SCRATCH ARRAY.
C
C JOB INTEGER.
C JOB CONTROLS THE COMPUTATION OF THE SINGULAR
C VECTORS. IT HAS THE DECIMAL EXPANSION AB
C WITH THE FOLLOWING MEANING
C
C A.EQ.0 DO NOT COMPUTE THE LEFT SINGULAR
C VECTORS.
C A.EQ.1 RETURN THE N LEFT SINGULAR VECTORS
C IN U.
C A.GE.2 RETURN THE FIRST MIN(N,P) SINGULAR
C VECTORS IN U.
C B.EQ.0 DO NOT COMPUTE THE RIGHT SINGULAR
C VECTORS.
C B.EQ.1 RETURN THE RIGHT SINGULAR VECTORS
C IN V.
C
C ON RETURN
C
C S DOUBLE PRECISION(MM), WHERE MM=MIN(N+1,P).
C THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
C SINGULAR VALUES OF X ARRANGED IN DESCENDING
C ORDER OF MAGNITUDE.
C
C E DOUBLE PRECISION(P)
C E ORDINARILY CONTAINS ZEROS. HOWEVER SEE THE
C DISCUSSION OF INFO FOR EXCEPTIONS.
C
C U DOUBLE PRECISION(LDU,K), WHERE LDU.GE.N. IF
C JOBA.EQ.1 THEN K.EQ.N, IF JOBA.GE.2
C THEN K.EQ.MIN(N,P).
C U CONTAINS THE MATRIX OF LEFT SINGULAR VECTORS.
C U IS NOT REFERENCED IF JOBA.EQ.0. IF N.LE.P
C OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
C IN THE SUBROUTINE CALL.
C
C V DOUBLE PRECISION(LDV,P), WHERE LDV.GE.P.
C V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C V IS NOT REFERENCED IF JOB.EQ.0. IF P.LE.N,
C THEN V MAY BE IDENTIFIED WITH X IN THE
C SUBROUTINE CALL.
C
C INFO INTEGER.
C THE SINGULAR VALUES (AND THEIR CORRESPONDING
C SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
C ARE CORRECT (HERE M=MIN(N,P)). THUS IF
C INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
C VECTORS ARE CORRECT. IN ANY EVENT, THE MATRIX
C B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
C WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
C ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
C IS THE TRANSPOSE OF U). THUS THE SINGULAR
C VALUES OF X AND B ARE THE SAME.
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C DSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
C
C EXTERNAL DROT
C BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG
C FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT
C
C INTERNAL VARIABLES
C
INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
* MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
DOUBLE PRECISION DDOT,T
DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,DNRM2,SCALE,SHIFT,SL,SM,SN,
* SMM1,T1,TEST,ZTEST
LOGICAL WANTU,WANTV
C
C
C SET THE MAXIMUM NUMBER OF ITERATIONS.
C
MAXIT = 30
C
C DETERMINE WHAT IS TO BE COMPUTED.
C
WANTU = .FALSE.
WANTV = .FALSE.
JOBU = MOD(JOB,100)/10
NCU = N
IF (JOBU .GT. 1) NCU = MIN0(N,P)
IF (JOBU .NE. 0) WANTU = .TRUE.
IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
C
C REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
C IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
C
INFO = 0
NCT = MIN0(N-1,P)
NRT = MAX0(0,MIN0(P-2,N))
LU = MAX0(NCT,NRT)
IF (LU .LT. 1) GO TO 170
DO 160 L = 1, LU
LP1 = L + 1
IF (L .GT. NCT) GO TO 20
C
C COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
C PLACE THE L-TH DIAGONAL IN S(L).
C
S(L) = DNRM2(N-L+1,X(L,L),1)
IF (S(L) .EQ. 0.0D0) GO TO 10
IF (X(L,L) .NE. 0.0D0) S(L) = DSIGN(S(L),X(L,L))
CALL DSCAL(N-L+1,1.0D0/S(L),X(L,L),1)
X(L,L) = 1.0D0 + X(L,L)
10 CONTINUE
S(L) = -S(L)
20 CONTINUE
IF (P .LT. LP1) GO TO 50
DO 40 J = LP1, P
IF (L .GT. NCT) GO TO 30
IF (S(L) .EQ. 0.0D0) GO TO 30
C
C APPLY THE TRANSFORMATION.
C
T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
30 CONTINUE
C
C PLACE THE L-TH ROW OF X INTO E FOR THE
C SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
C
E(J) = X(L,J)
40 CONTINUE
50 CONTINUE
IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
C
C PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
C MULTIPLICATION.
C
DO 60 I = L, N
U(I,L) = X(I,L)
60 CONTINUE
70 CONTINUE
IF (L .GT. NRT) GO TO 150
C
C COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
C L-TH SUPER-DIAGONAL IN E(L).
C
E(L) = DNRM2(P-L,E(LP1),1)
IF (E(L) .EQ. 0.0D0) GO TO 80
IF (E(LP1) .NE. 0.0D0) E(L) = DSIGN(E(L),E(LP1))
CALL DSCAL(P-L,1.0D0/E(L),E(LP1),1)
E(LP1) = 1.0D0 + E(LP1)
80 CONTINUE
E(L) = -E(L)
IF (LP1 .GT. N .OR. E(L) .EQ. 0.0D0) GO TO 120
C
C APPLY THE TRANSFORMATION.
C
DO 90 I = LP1, N
WORK(I) = 0.0D0
90 CONTINUE
DO 100 J = LP1, P
CALL DAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
100 CONTINUE
DO 110 J = LP1, P
CALL DAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
110 CONTINUE
120 CONTINUE
IF (.NOT.WANTV) GO TO 140
C
C PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
C BACK MULTIPLICATION.
C
DO 130 I = LP1, P
V(I,L) = E(I)
130 CONTINUE
140 CONTINUE
150 CONTINUE
160 CONTINUE
170 CONTINUE
C
C SET UP THE FINAL BIDIAGONAL MATRIX OF ORDER M.
C
M = MIN0(P,N+1)
NCTP1 = NCT + 1
NRTP1 = NRT + 1
IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
IF (N .LT. M) S(M) = 0.0D0
IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
E(M) = 0.0D0
C
C IF REQUIRED, GENERATE U.
C
IF (.NOT.WANTU) GO TO 300
IF (NCU .LT. NCTP1) GO TO 200
DO 190 J = NCTP1, NCU
DO 180 I = 1, N
U(I,J) = 0.0D0
180 CONTINUE
U(J,J) = 1.0D0
190 CONTINUE
200 CONTINUE
IF (NCT .LT. 1) GO TO 290
DO 280 LL = 1, NCT
L = NCT - LL + 1
IF (S(L) .EQ. 0.0D0) GO TO 250
LP1 = L + 1
IF (NCU .LT. LP1) GO TO 220
DO 210 J = LP1, NCU
T = -DDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
CALL DAXPY(N-L+1,T,U(L,L),1,U(L,J),1)
210 CONTINUE
220 CONTINUE
CALL DSCAL(N-L+1,-1.0D0,U(L,L),1)
U(L,L) = 1.0D0 + U(L,L)
LM1 = L - 1
IF (LM1 .LT. 1) GO TO 240
DO 230 I = 1, LM1
U(I,L) = 0.0D0
230 CONTINUE
240 CONTINUE
GO TO 270
250 CONTINUE
DO 260 I = 1, N
U(I,L) = 0.0D0
260 CONTINUE
U(L,L) = 1.0D0
270 CONTINUE
280 CONTINUE
290 CONTINUE
300 CONTINUE
C
C IF IT IS REQUIRED, GENERATE V.
C
IF (.NOT.WANTV) GO TO 350
DO 340 LL = 1, P
L = P - LL + 1
LP1 = L + 1
IF (L .GT. NRT) GO TO 320
IF (E(L) .EQ. 0.0D0) GO TO 320
DO 310 J = LP1, P
T = -DDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
CALL DAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1)
310 CONTINUE
320 CONTINUE
DO 330 I = 1, P
V(I,L) = 0.0D0
330 CONTINUE
V(L,L) = 1.0D0
340 CONTINUE
350 CONTINUE
C
C MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
C
MM = M
ITER = 0
360 CONTINUE
C
C QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
C
C ...EXIT
IF (M .EQ. 0) GO TO 620
C
C IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
C FLAG AND RETURN.
C
IF (ITER .LT. MAXIT) GO TO 370
INFO = M
C ......EXIT
GO TO 620
370 CONTINUE
C
C THIS SECTION OF THE PROGRAM INSPECTS FOR
C NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS. ON
C COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
C
C KASE = 1 IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
C KASE = 2 IF S(L) IS NEGLIGIBLE AND L.LT.M
C KASE = 3 IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
C S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
C KASE = 4 IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
C
DO 390 LL = 1, M
L = M - LL
C ...EXIT
IF (L .EQ. 0) GO TO 400
TEST = DABS(S(L)) + DABS(S(L+1))
ZTEST = TEST + DABS(E(L))
IF (ZTEST .NE. TEST) GO TO 380
E(L) = 0.0D0
C ......EXIT
GO TO 400
380 CONTINUE
390 CONTINUE
400 CONTINUE
IF (L .NE. M - 1) GO TO 410
KASE = 4
GO TO 480
410 CONTINUE
LP1 = L + 1
MP1 = M + 1
DO 430 LLS = LP1, MP1
LS = M - LLS + LP1
C ...EXIT
IF (LS .EQ. L) GO TO 440
TEST = 0.0D0
IF (LS .NE. M) TEST = TEST + DABS(E(LS))
IF (LS .NE. L + 1) TEST = TEST + DABS(E(LS-1))
ZTEST = TEST + DABS(S(LS))
IF (ZTEST .NE. TEST) GO TO 420
S(LS) = 0.0D0
C ......EXIT
GO TO 440
420 CONTINUE
430 CONTINUE
440 CONTINUE
IF (LS .NE. L) GO TO 450
KASE = 3
GO TO 470
450 CONTINUE
IF (LS .NE. M) GO TO 460
KASE = 1
GO TO 470
460 CONTINUE
KASE = 2
L = LS
470 CONTINUE
480 CONTINUE
L = L + 1
C
C PERFORM THE TASK INDICATED BY KASE.
C
GO TO (490,520,540,570), KASE
C
C DEFLATE NEGLIGIBLE S(M).
C
490 CONTINUE
MM1 = M - 1
F = E(M-1)
E(M-1) = 0.0D0
DO 510 KK = L, MM1
K = MM1 - KK + L
T1 = S(K)
CALL DROTG(T1,F,CS,SN)
S(K) = T1
IF (K .EQ. L) GO TO 500
F = -SN*E(K-1)
E(K-1) = CS*E(K-1)
500 CONTINUE
IF (WANTV) CALL DROT(P,V(1,K),1,V(1,M),1,CS,SN)
510 CONTINUE
GO TO 610
C
C SPLIT AT NEGLIGIBLE S(L).
C
520 CONTINUE
F = E(L-1)
E(L-1) = 0.0D0
DO 530 K = L, M
T1 = S(K)
CALL DROTG(T1,F,CS,SN)
S(K) = T1
F = -SN*E(K)
E(K) = CS*E(K)
IF (WANTU) CALL DROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
530 CONTINUE
GO TO 610
C
C PERFORM ONE QR STEP.
C
540 CONTINUE
C
C CALCULATE THE SHIFT.
C
SCALE = DMAX1(DABS(S(M)),DABS(S(M-1)),DABS(E(M-1)),
* DABS(S(L)),DABS(E(L)))
SM = S(M)/SCALE
SMM1 = S(M-1)/SCALE
EMM1 = E(M-1)/SCALE
SL = S(L)/SCALE
EL = E(L)/SCALE
B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0
C = (SM*EMM1)**2
SHIFT = 0.0D0
IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 550
SHIFT = DSQRT(B**2+C)
IF (B .LT. 0.0D0) SHIFT = -SHIFT
SHIFT = C/(B + SHIFT)
550 CONTINUE
F = (SL + SM)*(SL - SM) - SHIFT
G = SL*EL
C
C CHASE ZEROS.
C
MM1 = M - 1
DO 560 K = L, MM1
CALL DROTG(F,G,CS,SN)
IF (K .NE. L) E(K-1) = F
F = CS*S(K) + SN*E(K)
E(K) = CS*E(K) - SN*S(K)
G = SN*S(K+1)
S(K+1) = CS*S(K+1)
IF (WANTV) CALL DROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
CALL DROTG(F,G,CS,SN)
S(K) = F
F = CS*E(K) + SN*S(K+1)
S(K+1) = -SN*E(K) + CS*S(K+1)
G = SN*E(K+1)
E(K+1) = CS*E(K+1)
IF (WANTU .AND. K .LT. N)
* CALL DROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
560 CONTINUE
E(M-1) = F
ITER = ITER + 1
GO TO 610
C
C CONVERGENCE.
C
570 CONTINUE
C
C MAKE THE SINGULAR VALUE POSITIVE.
C
IF (S(L) .GE. 0.0D0) GO TO 580
S(L) = -S(L)
IF (WANTV) CALL DSCAL(P,-1.0D0,V(1,L),1)
580 CONTINUE
C
C ORDER THE SINGULAR VALUE.
C
590 IF (L .EQ. MM) GO TO 600
C ...EXIT
IF (S(L) .GE. S(L+1)) GO TO 600
T = S(L)
S(L) = S(L+1)
S(L+1) = T
IF (WANTV .AND. L .LT. P)
* CALL DSWAP(P,V(1,L),1,V(1,L+1),1)
IF (WANTU .AND. L .LT. N)
* CALL DSWAP(N,U(1,L),1,U(1,L+1),1)
L = L + 1
GO TO 590
600 CONTINUE
ITER = 0
M = M - 1
610 CONTINUE
GO TO 360
620 CONTINUE
RETURN
END
C
C ------------------------------------------------------------------
C
C TITLE: RS
C
C
C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
C OF A REAL SYMMETRIC MATRIX.
C
C ON INPUT:
C
C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT;
C
C N IS THE ORDER OF THE MATRIX A;
C
C A CONTAINS THE REAL SYMMETRIC MATRIX;
C
C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
C ONLY EIGENVALUES ARE DESIRED; OTHERWISE IT IS SET TO
C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
C
C ON OUTPUT:
C
C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER;
C
C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO;
C
C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN
C ERROR COMPLETION CODE DESCRIBED IN SECTION 2B OF THE
C DOCUMENTATION. THE NORMAL COMPLETION CODE IS ZERO;
C
C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR)
C
INTEGER N,NM,IERR,MATZ
DOUBLE PRECISION A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
C
IF (N .LE. NM) GO TO 10
IERR = 10 * N
GO TO 50
C
10 IF (MATZ .NE. 0) GO TO 20
C
C::::: FIND EIGENVALUES ONLY ::::::::::
CALL TRED1(NM,N,A,W,FV1,FV2)
CALL TQLRAT(N,W,FV2,IERR)
GO TO 50
C
C::::: FIND BOTH EIGENVALUES AND EIGENVECTORS
C:::::
20 CALL TRED2(NM,N,A,W,FV1,Z)
CALL TQL2(NM,N,W,FV1,Z,IERR)
50 RETURN
C
C::::: LAST CARD OF RS
C::
END
C
C ------------------------------------------------------------------
C
SUBROUTINE TRED1(NM,N,A,D,E,E2)
C
INTEGER I,J,K,L,N,II,NM,JP1
DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N)
DOUBLE PRECISION F,G,H,SCALE
DOUBLE PRECISION DSQRT,DABS,DSIGN
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1,
C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX
C TO A SYMMETRIC TRIDIAGONAL MATRIX USING
C ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C ON INPUT:
C
C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT;
C
C N IS THE ORDER OF THE MATRIX;
C
C A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE
C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C ON OUTPUT:
C
C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
C FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
C TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED;
C
C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX;
C
C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO;
C
C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
DO 100 I = 1, N
100 D(I) = A(I,I)
C
C::::: FOR I=N STEP -1 UNTIL 1 DO --
C::::::
DO 300 II = 1, N
I = N + 1 - II
L = I - 1
H = 0.0D0
SCALE = 0.0D0
IF (L .LT. 1) GO TO 130
C
C::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED)
C::::
DO 120 K = 1, L
120 SCALE = SCALE + DABS(A(I,K))
C
IF (SCALE .NE. 0.0D0) GO TO 140
130 E(I) = 0.0D0
E2(I) = 0.0D0
GO TO 290
C
140 DO 150 K = 1, L
A(I,K) = A(I,K) / SCALE
H = H + A(I,K) * A(I,K)
150 CONTINUE
C
E2(I) = SCALE * SCALE * H
F = A(I,L)
G = -DSIGN(DSQRT(H),F)
E(I) = SCALE * G
H = H - F * G
A(I,L) = F - G
IF (L .EQ. 1) GO TO 270
F = 0.0D0
C
DO 240 J = 1, L
G = 0.0D0
C
C::::: FORM ELEMENT OF A*U
C::::::
DO 180 K = 1, J
180 G = G + A(J,K) * A(I,K)
C
JP1 = J + 1
IF (L .LT. JP1) GO TO 220
C
DO 200 K = JP1, L
200 G = G + A(K,J) * A(I,K)
C
C::::: FORM ELEMENT OF P
C::::
220 E(J) = G / H
F = F + E(J) * A(I,J)
240 CONTINUE
C
H = F / (H + H)
C
C::::: FORM REDUCED A
C:
DO 260 J = 1, L
F = A(I,J)
G = E(J) - H * F
E(J) = G
C
DO 260 K = 1, J
A(J,K) = A(J,K) - F * E(K) - G * A(I,K)
260 CONTINUE
C
270 DO 280 K = 1, L
280 A(I,K) = SCALE * A(I,K)
C
290 H = D(I)
D(I) = A(I,I)
A(I,I) = H
300 CONTINUE
C
RETURN
C
C::::: LAST CARD OF TRED1
C:::::
END
C
C ------------------------------------------------------------------
C
SUBROUTINE TQLRAT(N,D,E2,IERR)
C
INTEGER I,J,L,M,N,II,L1,MML,IERR
DOUBLE PRECISION D(N),E2(N)
DOUBLE PRECISION B,C,F,G,H,P,R,S,T,MACHEP
DOUBLE PRECISION DSQRT,DABS,DSIGN
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
C ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
C
C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
C TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD.
C
C ON INPUT:
C
C N IS THE ORDER OF THE MATRIX;
C
C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX;
C
C E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE
C INPUT MATRIX IN ITS LAST N-1 POSITIONS. E2(1) IS ARBITRARY.
C
C ON OUTPUT:
C
C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
C THE SMALLEST EIGENVALUES;
C
C E2 HAS BEEN DESTROYED;
C
C IERR IS SET TO
C ZERO FOR NORMAL RETURN,
C J IF THE J-TH EIGENVALUE HAS NOT BEEN
C DETERMINED AFTER 30 ITERATIONS.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
C
C::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C AN APPROXIMATION TO MACHEP IS COMPUTED BELOW.
C
C::::::
MACHEP = 1.D0 / 2.D0
DO 2 I=1,500
T = 1.D0 + MACHEP
IF (T .EQ. 1.D0) GO TO 3
MACHEP = MACHEP / 2.D0
2 CONTINUE
C
3 IERR = 0
IF (N .EQ. 1) GO TO 1001
C
DO 100 I = 2, N
100 E2(I-1) = E2(I)
C
F = 0.0D0
B = 0.0D0
E2(N) = 0.0D0
C
DO 290 L = 1, N
J = 0
H = MACHEP * (DABS(D(L)) + DSQRT(E2(L)))
IF (B .GT. H) GO TO 105
B = H
C = B * B
C
C::::: LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT
C
105 DO 110 M = L, N
IF (E2(M) .LE. C) GO TO 120
C
C::::: E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C THROUGH THE BOTTOM OF THE LOOP
C:::::::
110 CONTINUE
C
120 IF (M .EQ. L) GO TO 210
130 IF (J .EQ. 30) GO TO 1000
J = J + 1
C
C::::: FORM SHIFT
C:::::::
L1 = L + 1
S = DSQRT(E2(L))
G = D(L)
P = (D(L1) - G) / (2.0D0 * S)
R = DSQRT(P*P+1.0D0)
D(L) = S / (P + DSIGN(R,P))
H = G - D(L)
C
DO 140 I = L1, N
140 D(I) = D(I) - H
C
F = F + H
C
C::::: RATIONAL QL TRANSFORMATION
C:::
G = D(M)
IF (G .EQ. 0.0D0) G = B
H = G
S = 0.0D0
MML = M - L
C
C::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
DO 200 II = 1, MML
I = M - II
P = G * H
R = P + E2(I)
E2(I+1) = S * R
S = E2(I) / R
D(I+1) = H + S * (H + D(I))
G = D(I) - E2(I) / G
IF (G .EQ. 0.0D0) G = B
H = G * P / R
200 CONTINUE
C
E2(L) = S * G
D(L) = H
C
C::::: GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST
C
IF (H .EQ. 0.0D0) GO TO 210
IF (DABS(E2(L)) .LE. DABS(C/H)) GO TO 210
E2(L) = H * E2(L)
IF (E2(L) .NE. 0.0D0) GO TO 130
210 P = D(L) + F
C
C::::: ORDER EIGENVALUES
C::::
IF (L .EQ. 1) GO TO 250
C
C::::: FOR I=L STEP -1 UNTIL 2 DO --
C::::::
DO 230 II = 2, L
I = L + 2 - II
IF (P .GE. D(I-1)) GO TO 270
D(I) = D(I-1)
230 CONTINUE
C
250 I = 1
270 D(I) = P
290 CONTINUE
C
GO TO 1001
C
C::::: SET ERROR -- NO CONVERGENCE TO AN
C EIGENVALUE AFTER 30 ITERATIONS
C:::::::
1000 IERR = L
1001 RETURN
C
C::::: LAST CARD OF TQLRAT
C::::::
END
C
C ------------------------------------------------------------------
C
SUBROUTINE TRED2(NM,N,A,D,E,Z)
C
INTEGER I,J,K,L,N,II,NM,JP1
DOUBLE PRECISION A(NM,N),D(N),E(N),Z(NM,N)
DOUBLE PRECISION F,G,H,HH,SCALE
DOUBLE PRECISION DSQRT,DABS,DSIGN
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A
C SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
C ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C ON INPUT:
C
C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT;
C
C N IS THE ORDER OF THE MATRIX;
C
C A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE
C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C ON OUTPUT:
C
C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX;
C
C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO;
C
C Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
C PRODUCED IN THE REDUCTION;
C
C A AND Z MAY COINCIDE. IF DISTINCT, A IS UNALTERED.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
DO 100 I = 1, N
C
DO 100 J = 1, I
Z(I,J) = A(I,J)
100 CONTINUE
C
IF (N .EQ. 1) GO TO 320
C
C::::: FOR I=N STEP -1 UNTIL 2 DO --
C::::::
DO 300 II = 2, N
I = N + 2 - II
L = I - 1
H = 0.0D0
SCALE = 0.0D0
IF (L .LT. 2) GO TO 130
C
C::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED)
C::::
DO 120 K = 1, L
120 SCALE = SCALE + DABS(Z(I,K))
C
IF (SCALE .NE. 0.0D0) GO TO 140
130 E(I) = Z(I,L)
GO TO 290
C
140 DO 150 K = 1, L
Z(I,K) = Z(I,K) / SCALE
H = H + Z(I,K) * Z(I,K)
150 CONTINUE
C
F = Z(I,L)
G = -DSIGN(DSQRT(H),F)
E(I) = SCALE * G
H = H - F * G
Z(I,L) = F - G
F = 0.0D0
C
DO 240 J = 1, L
Z(J,I) = Z(I,J) / H
G = 0.0D0
C
C::::: FORM ELEMENT OF A*U
C::::::
DO 180 K = 1, J
180 G = G + Z(J,K) * Z(I,K)
C
JP1 = J + 1
IF (L .LT. JP1) GO TO 220
C
DO 200 K = JP1, L
200 G = G + Z(K,J) * Z(I,K)
C
C::::: FORM ELEMENT OF P
C::::
220 E(J) = G / H
F = F + E(J) * Z(I,J)
240 CONTINUE
C
HH = F / (H + H)
C
C::::: FORM REDUCED A
C:
DO 260 J = 1, L
F = Z(I,J)
G = E(J) - HH * F
E(J) = G
C
DO 260 K = 1, J
Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
260 CONTINUE
C
290 D(I) = H
300 CONTINUE
C
320 D(1) = 0.0D0
E(1) = 0.0D0
C
C::::: ACCUMULATION OF TRANSFORMATION MATRICES
C::::::
DO 500 I = 1, N
L = I - 1
IF (D(I) .EQ. 0.0D0) GO TO 380
C
DO 360 J = 1, L
G = 0.0D0
C
DO 340 K = 1, L
340 G = G + Z(I,K) * Z(K,J)
C
DO 360 K = 1, L
Z(K,J) = Z(K,J) - G * Z(K,I)
360 CONTINUE
C
380 D(I) = Z(I,I)
Z(I,I) = 1.0D0
IF (L .LT. 1) GO TO 500
C
DO 400 J = 1, L
Z(I,J) = 0.0D0
Z(J,I) = 0.0D0
400 CONTINUE
C
500 CONTINUE
C
RETURN
C
C::::: LAST CARD OF TRED2
C:::::
END
C
C ------------------------------------------------------------------
C
SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
C
INTEGER I,J,K,L,M,N,II,L1,NM,MML,IERR
DOUBLE PRECISION D(N),E(N),Z(NM,N)
DOUBLE PRECISION B,C,F,G,H,P,R,S,T,MACHEP
DOUBLE PRECISION DSQRT,DABS,DSIGN
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
C NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
C WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
C
C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS
C FULL MATRIX TO TRIDIAGONAL FORM.
C
C ON INPUT:
C
C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C DIMENSION STATEMENT;
C
C N IS THE ORDER OF THE MATRIX;
C
C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX;
C
C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY;
C
C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS
C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
C THE IDENTITY MATRIX.
C
C ON OUTPUT:
C
C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
C UNORDERED FOR INDICES 1,2,...,IERR-1;
C
C E HAS BEEN DESTROYED;
C
C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE,
C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
C EIGENVALUES;
C
C IERR IS SET TO
C ZERO FOR NORMAL RETURN,
C J IF THE J-TH EIGENVALUE HAS NOT BEEN
C DETERMINED AFTER 30 ITERATIONS.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C ------------------------------------------------------------------
C
C
C::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C AN APPROXIMATION TO MACHEP IS COMPUTED BELOW.
C
C::::::
MACHEP = 1.D0 / 2.D0
DO 2 I=1,500
T = 1.D0 + MACHEP
IF (T .EQ. 1.D0) GO TO 3
MACHEP = MACHEP / 2.D0
2 CONTINUE
C
3 IERR = 0
IF (N .EQ. 1) GO TO 1001
C
DO 100 I = 2, N
100 E(I-1) = E(I)
C
F = 0.0D0
B = 0.0D0
E(N) = 0.0D0
C
DO 240 L = 1, N
J = 0
H = MACHEP * (DABS(D(L)) + DABS(E(L)))
IF (B .LT. H) B = H
C
C::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT
C::
DO 110 M = L, N
IF (DABS(E(M)) .LE. B) GO TO 120
C
C::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C THROUGH THE BOTTOM OF THE LOOP
C:::::::
110 CONTINUE
C
120 IF (M .EQ. L) GO TO 220
130 IF (J .EQ. 30) GO TO 1000
J = J + 1
C
C::::: FORM SHIFT
C:::::::
L1 = L + 1
G = D(L)
P = (D(L1) - G) / (2.0D0 * E(L))
R = DSQRT(P*P+1.0D0)
D(L) = E(L) / (P + DSIGN(R,P))
H = G - D(L)
C
DO 140 I = L1, N
140 D(I) = D(I) - H
C
F = F + H
C
C::::: QL TRANSFORMATION
C::::
P = D(M)
C = 1.0D0
S = 0.0D0
MML = M - L
C
C::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
DO 200 II = 1, MML
I = M - II
G = C * E(I)
H = C * P
IF (DABS(P) .LT. DABS(E(I))) GO TO 150
C = E(I) / P
R = DSQRT(C*C+1.0D0)
E(I+1) = S * P * R
S = C / R
C = 1.0D0 / R
GO TO 160
150 C = P / E(I)
R = DSQRT(C*C+1.0D0)
E(I+1) = S * E(I) * R
S = 1.0D0 / R
C = C * S
160 P = C * D(I) - S * G
D(I+1) = H + S * (C * G + S * D(I))
C
C::::: FORM VECTOR ::::::::::
DO 180 K = 1, N
H = Z(K,I+1)
Z(K,I+1) = S * Z(K,I) + C * H
Z(K,I) = C * Z(K,I) - S * H
180 CONTINUE
C
200 CONTINUE
C
E(L) = S * P
D(L) = C * P
IF (DABS(E(L)) .GT. B) GO TO 130
220 D(L) = D(L) + F
240 CONTINUE
C
C::::: ORDER EIGENVALUES AND EIGENVECTORS
C:
DO 300 II = 2, N
I = II - 1
K = I
P = D(I)
C
DO 260 J = II, N
IF (D(J) .GE. P) GO TO 260
K = J
P = D(J)
260 CONTINUE
C
IF (K .EQ. I) GO TO 300
D(K) = D(I)
D(I) = P
C
DO 280 J = 1, N
P = Z(J,I)
Z(J,I) = Z(J,K)
Z(J,K) = P
280 CONTINUE
C
300 CONTINUE
C
GO TO 1001
C
C::::: SET ERROR -- NO CONVERGENCE TO AN
C EIGENVALUE AFTER 30 ITERATIONS
C:::::::
1000 IERR = L
1001 RETURN
C
C::::: LAST CARD OF TQL2
C::::
END
SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
C
C CONSTANT TIMES A VECTOR PLUS A VECTOR.
C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C JACK DONGARRA, LINPACK, 3/11/78.
C
DOUBLE PRECISION DX(1),DY(1),DA
INTEGER I,INCX,INCY,IXIY,M,MP1,N
C
IF(N.LE.0)RETURN
IF (DA .EQ. 0.0D0) RETURN
IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C NOT EQUAL TO 1
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DY(IY) = DY(IY) + DA*DX(IX)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C CLEAN-UP LOOP
C
20 M = MOD(N,4)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
DY(I) = DY(I) + DA*DX(I)
30 CONTINUE
IF( N .LT. 4 ) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,4
DY(I) = DY(I) + DA*DX(I)
DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
50 CONTINUE
RETURN
END
SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
C
C COPIES A VECTOR, X, TO A VECTOR, Y.
C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C JACK DONGARRA, LINPACK, 3/11/78.
C
DOUBLE PRECISION DX(1),DY(1)
INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
IF(N.LE.0)RETURN
IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C NOT EQUAL TO 1
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DY(IY) = DX(IX)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C CLEAN-UP LOOP
C
20 M = MOD(N,7)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
DY(I) = DX(I)
30 CONTINUE
IF( N .LT. 7 ) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,7
DY(I) = DX(I)
DY(I + 1) = DX(I + 1)
DY(I + 2) = DX(I + 2)
DY(I + 3) = DX(I + 3)
DY(I + 4) = DX(I + 4)
DY(I + 5) = DX(I + 5)
DY(I + 6) = DX(I + 6)
50 CONTINUE
RETURN
END
DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
C
C FORMS THE DOT PRODUCT OF TWO VECTORS.
C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C JACK DONGARRA, LINPACK, 3/11/78.
C
DOUBLE PRECISION DX(1),DY(1),DTEMP
INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
DDOT = 0.0D0
DTEMP = 0.0D0
IF(N.LE.0)RETURN
IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C NOT EQUAL TO 1
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DTEMP = DTEMP + DX(IX)*DY(IY)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
DDOT = DTEMP
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C CLEAN-UP LOOP
C
20 M = MOD(N,5)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
DTEMP = DTEMP + DX(I)*DY(I)
30 CONTINUE
IF( N .LT. 5 ) GO TO 60
40 MP1 = M + 1
DO 50 I = MP1,N,5
DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
* DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
50 CONTINUE
60 DDOT = DTEMP
RETURN
END
DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
INTEGER NEXT
DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
DATA ZERO, ONE /0.0D0, 1.0D0/
C
C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
C INCREMENT INCX .
C IF N .LE. 0 RETURN WITH RESULT = 0.
C IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C C.L.LAWSON, 1978 JAN 08
C
C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE
C HOPEFULLY APPLICABLE TO ALL MACHINES.
C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES.
C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES.
C WHERE
C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT)
C V = LARGEST NO. (OVERFLOW LIMIT)
C
C BRIEF OUTLINE OF ALGORITHM..
C
C PHASE 1 SCANS ZERO COMPONENTS.
C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C VALUES FOR CUTLO AND CUTHI..
C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE
C UNIVAC AND DEC AT 2**(-103)
C THUS CUTLO = 2**(-51) = 4.44089E-16
C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C THUS CUTHI = 2**(63.5) = 1.30438E19
C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C THUS CUTLO = 2**(-33.5) = 8.23181D-11
C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19
C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
C
IF(N .GT. 0) GO TO 10
DNRM2 = ZERO
GO TO 300
C
10 ASSIGN 30 TO NEXT
SUM = ZERO
NN = N * INCX
C BEGIN MAIN LOOP
I = 1
20 GO TO NEXT,(30, 50, 70, 110)
30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
ASSIGN 50 TO NEXT
XMAX = ZERO
C
C PHASE 1. SUM IS ZERO
C
50 IF( DX(I) .EQ. ZERO) GO TO 200
IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
C
C PREPARE FOR PHASE 2.
ASSIGN 70 TO NEXT
GO TO 105
C
C PREPARE FOR PHASE 4.
C
100 I = J
ASSIGN 110 TO NEXT
SUM = (SUM / DX(I)) / DX(I)
105 XMAX = DABS(DX(I))
GO TO 115
C
C PHASE 2. SUM IS SMALL.
C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
C
C COMMON CODE FOR PHASES 2 AND 4.
C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
C
110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
SUM = ONE + SUM * (XMAX / DX(I))**2
XMAX = DABS(DX(I))
GO TO 200
C
115 SUM = SUM + (DX(I)/XMAX)**2
GO TO 200
C
C
C PREPARE FOR PHASE 3.
C
75 SUM = (SUM * XMAX) * XMAX
C
C
C FOR REAL OR D.P. SET HITEST = CUTHI/N
C FOR COMPLEX SET HITEST = CUTHI/(2*N)
C
85 HITEST = CUTHI/FLOAT( N )
C
C PHASE 3. SUM IS MID-RANGE. NO SCALING.
C
DO 95 J =I,NN,INCX
IF(DABS(DX(J)) .GE. HITEST) GO TO 100
95 SUM = SUM + DX(J)**2
DNRM2 = DSQRT( SUM )
GO TO 300
C
200 CONTINUE
I = I + INCX
IF ( I .LE. NN ) GO TO 20
C
C END OF MAIN LOOP.
C
C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
DNRM2 = XMAX * DSQRT(SUM)
300 CONTINUE
RETURN
END
SUBROUTINE DROT (N,DX,INCX,DY,INCY,C,S)
00000010
C
C APPLIES A PLANE ROTATION.
C JACK DONGARRA, LINPACK, 3/11/78.
C
DOUBLE PRECISION DX(1),DY(1),DTEMP,C,S
INTEGER I,INCX,INCY,IX,IY,N
C
IF(N.LE.0)RETURN
IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C TO 1
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DTEMP = C*DX(IX) + S*DY(IY)
DY(IY) = C*DY(IY) - S*DX(IX)
DX(IX) = DTEMP
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1
C
20 DO 30 I = 1,N
DTEMP = C*DX(I) + S*DY(I)
DY(I) = C*DY(I) - S*DX(I)
DX(I) = DTEMP
30 CONTINUE
RETURN
END
SUBROUTINE DROTG(DA,DB,C,S)
00000010
C
C CONSTRUCT GIVENS PLANE ROTATION.
C JACK DONGARRA, LINPACK, 3/11/78.
C
DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z
C
ROE = DB
IF( DABS(DA) .GT. DABS(DB) ) ROE = DA
SCALE = DABS(DA) + DABS(DB)
IF( SCALE .NE. 0.0D0 ) GO TO 10
C = 1.0D0
S = 0.0D0
R = 0.0D0
GO TO 20
10 R = SCALE*DSQRT((DA/SCALE)**2 + (DB/SCALE)**2)
R = DSIGN(1.0D0,ROE)*R
C = DA/R
S = DB/R
20 Z = 1.0D0
IF( DABS(DA) .GT. DABS(DB) ) Z = S
IF( DABS(DB) .GE. DABS(DA) .AND. C .NE. 0.0D0 ) Z = 1.0D0/C
DA = R
DB = Z
RETURN
END
SUBROUTINE DSCAL(N,DA,DX,INCX)
C
C SCALES A VECTOR BY A CONSTANT.
C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
C JACK DONGARRA, LINPACK, 3/11/78.
C
DOUBLE PRECISION DA,DX(1)
INTEGER I,INCX,M,MP1,N,NINCX
C
IF(N.LE.0)RETURN
IF(INCX.EQ.1)GO TO 20
C
C CODE FOR INCREMENT NOT EQUAL TO 1
C
NINCX = N*INCX
DO 10 I = 1,NINCX,INCX
DX(I) = DA*DX(I)
10 CONTINUE
RETURN
C
C CODE FOR INCREMENT EQUAL TO 1
C
C
C CLEAN-UP LOOP
C
20 M = MOD(N,5)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
DX(I) = DA*DX(I)
30 CONTINUE
IF( N .LT. 5 ) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,5
DX(I) = DA*DX(I)
DX(I + 1) = DA*DX(I + 1)
DX(I + 2) = DA*DX(I + 2)
DX(I + 3) = DA*DX(I + 3)
DX(I + 4) = DA*DX(I + 4)
50 CONTINUE
RETURN
END
SUBROUTINE DSWAP (N,DX,INCX,DY,INCY)
C
C INTERCHANGES TWO VECTORS.
C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
C JACK DONGARRA, LINPACK, 3/11/78.
C
DOUBLE PRECISION DX(1),DY(1),DTEMP
INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
IF(N.LE.0)RETURN
IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C TO 1
C
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DTEMP = DX(IX)
DX(IX) = DY(IY)
DY(IY) = DTEMP
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C CLEAN-UP LOOP
C
20 M = MOD(N,3)
IF( M .EQ. 0 ) GO TO 40
DO 30 I = 1,M
DTEMP = DX(I)
DX(I) = DY(I)
DY(I) = DTEMP
30 CONTINUE
IF( N .LT. 3 ) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,3
DTEMP = DX(I)
DX(I) = DY(I)
DY(I) = DTEMP
DTEMP = DX(I + 1)
DX(I + 1) = DY(I + 1)
DY(I + 1) = DTEMP
DTEMP = DX(I + 2)
DX(I + 2) = DY(I + 2)
DY(I + 2) = DTEMP
50 CONTINUE
RETURN
END
LONGLEY DATA EXAMPLE FROM JASA 1967 62 819-841
16 6 10 30 1 2
83.0 234289 2356 1590 107608 1947 60323
88.5 259426 2325 1456 108632 1948 61122
88.2 258054 3682 1616 109773 1949 60171
89.5 284599 3351 1650 110929 1950 61187
96.2 328975 2099 3099 112075 1951 63221
98.1 346999 1932 3594 113270 1952 63639
99.0 365385 1870 3547 115094 1953 64989
100.0 363112 3578 3350 116219 1954 63761
101.2 397469 2904 3048 117388 1955 66019
104.6 419180 2822 2857 118734 1956 67857
108.4 442769 2936 2798 120445 1957 68169
110.8 444546 4681 2637 121950 1958 66513
112.6 482704 3813 2552 123366 1959 68655
114.2 502601 3931 2514 125368 1960 69564
115.7 518173 4806 2572 127852 1961 69331
116.9 554894 4007 2827 130081 1962 70551
16 6 21 0 3 2